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PREFACE 


This brochure contains some 30 programs written for a 
specific programmable pocket calculator, the HP-25. 
These programs implement algorithms in number theory, 
equation solving, algebraic stability theory, numeri- 
cal integration, as well as for the evaluation of 
special functions such as the gamma function, various 
Bessel functions, and the Riemann zeta function. By 
means of the flow diagrams and the detailed descrip- 
tions that are provided, the programs are easily 
adapted to run on any calculator of comparable capa- 
CLty. 

The purpose to be served by these programs is di- 
dactic as well as practical. 

The immediate didactic purpose is to enable stu- 
dents of numerical and computational analysis to gain 
first-hand experience with some modern techniques of 
scientific computing. I have used a number of the 
programs in the classroom to demonstrate the numeri- 
cal performance of such techniques, using data provi- 
ded spontaneously by the students. Such demonstra- 
tions greatly increase the practical, "know-how" con- 
tent of numerical analysis courses. They also elimi- 
nate the suspicion that the instructor uses examples 
that are rigged to show an algorithm in an especially 
favorable light. 

My purpose is didactic in another, potentially 
even more important sense. Programmable calculators 
have now been in use for almost a third of a century. 
Yet many members of the scientific community, al- 
though they use mathematics in their work, have re- 
mained ignorant of the essentials of modern scienti- 
fic computation. This even holds for some mathemati- 
Cians, if they did not choose to specialize in nume- 
rical analysis or computer science. All this is 
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destined to change with the advent of the program- 
mable pocket calculator. By providing the basic fea- 
tures of branching and iteration, these amazing in- 
struments put automatic computation within the reach, 
in the privacy of one's own study, of anybody who re- 
members his basic calculus. No special programming 
language must be learned to operate one of these de- 
vices; to its owner, the access is immediate, and 
there are no waiting times beyond those required by 
the computation itself. An error in the program or in 
the data can be corrected immediately. Thus these 
calculators are interactive to a degree that up to 
now was enjoyed only by the privileged few who were 
connected to an expensive computer by their own ter- 
minal. By publishing my programs I hope to bring some 
of the flavor of modern automatic computation to 
those who so far have been untouched by it. 

Finally, some of my programs, especially those 
dealing with special functions, also serve an emi- 
nently practical purpose. If I quickly need some va- 
lues of a higher transcendental function such as a 
Bessel function of nonintegral order, I know of no 
way to get them with less fuss than by using my 
little calculator. Similar programs could (and un- 
doubtedly will) be written for functions other than 
those represented in my collection; my programs mere- 
ly provide a small sample of what can be done in this 
area. 

Coming back to the didactic aspect, it must be 
admitted that not all areas of modern numerical com- 
putation can be demonstrated satisfactorily ona 
pocket calculator. Because of the limited memory ca- 
pacity, problems dealing with large sets of data, 
such as computations in linear algebra, optimization, 
or partial differential equations, cannot be solved 
on a pocket calculator as they would be solved ona 
large computer. Even in problems requiring small sets 
of data, the lack of memory space imposes limitations. 
Thus it happens, for instance, that my programs for 
determining the zeros or the stability properties of 
polynomials in most cases can deal only with polyno- 
mials of degree not exceeding four. Such programs ne- 
vertheless can serve as models - or "pilot programs" 
- for similar programs designed to deal with larger 
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sets of data on correspondingly larger computers. In- 
deed, some of my programs for the evaluation of spe- 
cial functions would be written in much the same way 
also for the very largest computers; the only diffe- 
rence here lies in the speed of execution. 

Naturally, this collection is not meant to be ex- 
haustive. The creative student can find ample room 
for the design of new programs, particularly in the 
areas of equation solving, numerical integration, and 
the evaluation of special functions. Furthermore, it 
is likely that many of the programs presented here 
could be improved in some way, say, by making them 
faster, easier to use, or more automatic. I will be 
genuinely pleased to learn about such improvements. 

Already when compiling these programs I have be- 
nefitted from the advice and help of numerous colle- 
gues, students, and friends. J. Waldvogel and D. D. 
Warner contributed ideas. W. Seewald and A. Stahli 
Suggested improvements in existing programs. 

E. Specker offered a solution to the problem of wri- 
ting viable power series programs. My wife, Marie- 
Louise Henrici, did most of the checking and proof- 
reading and eliminated numerous errors. The reader 
may judge for himself the excellent work of Brigitte 
Knecht who transformed my manuscript into camera- 
ready copy. The staff of John Wiley & Sons disposed 
of all problems of editing and production with pro- 
fessional know-how. To all these individuals, and to 
any other helpers who may not have been named here, 
I offer my heartfelt thanks. 


Peter Henrici 
Zurich, Switzerland 
February 1977 


CONTENTS 


INTRODUCTION 


Part 1. NUMBER THEORY 


Prime Factor Decomposition 
Euclidean Algorithm 
Rational Binomial Coefficients 


Continued Fraction Representation 
of Real Numbers 


Exact Continued Fractions for 
Quadratic Irrationalities 
Part 2. ITERATION 


Iteration 
Iteration with Aitken Acceleration 
Aitken-Steffensen Iteration 


Newton Iteration for Complex Roots 


Part 3. POLYNOMIALS 


Horner Algorithm 


67 


1x 


Xx Contents 


Newton's Method for Polynomials 73 
Bernoulli's Method for Single Dominant Zero 80 
Bernoulli's Method for Complex Conjugate 87 
Dominant Zeros 
Quotient-Difference Algorithm 95 
Routh Algorithm 101 
Schur-Cohn Algorithm I 107 
Schur-Cohn Algorithm II 114 
Part 4. POWER SERIES 121 
Reciprocal Power Series 123 
Power of Power Series 134 
Exponentiation of Power Series 144 
Part 5. INTEGRATION 153 
Numerical Integration with Step Refinement 155 
Romberg Algorithm 162 
Plana Summation Formula 170 
Differential Equation of First Order, 181 
Trapezoidal Method 
Autonomous Differential Equation of Second 189 
Order, First Derivative Absent 
Linear Second Order Differential Equation 196 
Part 6. SPECIAL FUNCTIONS 205 


Log-Arcsine Algorithm 207 


Contents 


Gamma Function 

Incomplete Gamma Function 

Error Function 

Complete Elliptic Integrals 

Bessel Functions, Integer Order 
Bessel Functions, Arbitrary Order 
Bessel Functions: Asymptotic Series 


Riemann Zeta Function on Critical Line 


INDEX 


xi 


214 
223 
233 
241 
247 
253 
260 
268 


277 


Computational Analysis 
with the 
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INTRODUCTION 
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This is not a textbook on the programming of pocket 
calculators. To learn how to operate your calculator, 
consult your manual. As a first introduction to the 
manifold uses of the pocket calculator, the recent 
book by Jon M. Smith, "Scientific analysis on the 
pocket calculator" (Wiley, 1975) can be recommended. 
Each program description is arranged according 


to the following scheme: 


. Purpose, 


. Method, 


1 

2 

3. Flow diagram, 
4, Storage and program, 

5. Operating instructions, 
6 


. Examples and timing. 


The statement of purpose usually is very brief. 

It is intended to state, in language as nontechnical 
as possible, what the program is good for. 

The description of the method is strictly confined 
to what is being done. For the analytical justifica- 
tion and for technical details, the reader is referred 
to the literature. References to books by the author 


are given in the following abbreviated form: 


4 Introduction 


ENA = "Elements of numerical analysis", 


Wiley, New York, 1964, 


ACCA I = "Applied and computational complex ana- 
lysis", Vol. I, Wiley, New York, 1974, 
ACCA II = "Applied and computational complex ana- 


lysis", Vol. II, Wiley, New York, 1977. 


The flow diagrams do not follow any particular 
standard format. The intent is simply to make the 
structure of a program visible at a glance. In some 
cases explanations are given for the technical de- 
tails of a program or for abbreviations that are 
used in the flow diagram. 

The section entitled storage and program, natural- 
ly, is the core part of each program description. We 


use the symbol R, to denote the storage register num- 


k 
bered k. Quantities that must be stored in the storage 


registers ahead of the computation are encased, thus: 


The program listing uses the key symbols as they are 
found on the keys of your HP-25 calculator, with the 
exception of multiplication, which is indicated by * 


x". The 


in order to avoid confusion with the letter 
listing is as compact as possible. For yellow and 
blue instructions it is understood that the keys "f" 
and "g" are to be pressed in advance; no ambiguity 


can arise by our abbreviated notation. The arrows > 
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in front of certain instruction numbers are not part 
of the program. They simply indicate an entrance from 
a GTO instruction. For those readers who would have 
liked to see more complete descriptions of the tech- 
nical details of a program, it must be said that such 
paraphrases can make for extremely tedious reading. 
The best way to understand the details of a program 
is for the user to try to write his own. He then will 
often realize why the particular limitations of the 
pocket calculator forced the author to deviate from 
the more straightforward ways in which such programs 
would be written for a larger computer. 

An attempt has been made to keep the operating 
instructions self-contained, except in a few obvious 
cases. Here we also give indications and explanations 
of possible failures of some programs. 

The examples are not intended to be particularly 
sophisticated. The user should try one of the simpler 
examples to see whether he has correctly loaded the 
program and followed the instructions. In some cases 
the examples indicate how the program would be used 
in the classroom, illustrating particular points of 


analysis. All timings given are approximate. 


Part 1 
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PRIME FACTOR DECOMPOSITION 


1. Purpose 


To produce, for any positive integer m < 10" Lis ‘de= 
composition in prime factors. Factors are to appear in 


increasing order. 


2. Method 


We divide m successively by the numbers 


that is, by 2, 3 and by all odd integers > 3 that are 
not divisible by 3, until either (a) d > Ym or (b) the 
remainder is zero. In case (a), m is prime; the calcu- 
lator shows m and then zero to indicate that the de- 
composition is terminated. In case (b), dis a divisor 
of m. The calculator shows d, replaces m by m/d, and 
starts again, using the last das its first divisor. 
To compute the d, we start with d = 2 and calcu- 


late the successor d' of d by the formula 
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V4 20 Ae a 


4 
G+ 3.4 € 5, 2f:d > 4 


Here the value of e€ alternates between +1 and -l, be- 


ginning with -l. 
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3. Flow Diagram 


12 


4. Storage and Program 


R 


0 


Number Theory 


~ 
A 
Nm OC F KR FP NN BF 
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oe Operating Instructions 


After loading the program, move the operating switch 


to RUN. Then press 
FIX O 


to get the nine-digit display of all integers. If 
prime factors of m are desired, load m into the X 


register. When 


PRGM 
R/S 


is pressed, the calculation will start, then stop by 


displaying smallest prime factor of m. When 
R/S 


is pressed, the next prime factor will be displayed, 
and so forth, until 0 is shown, indicating that all 


prime factors have been found. 


6. Examples and Timing 


[1] 36 = 2 * 2 * 3 * 3. Time required about 8 sec. 
71489 = 11 * 67 * 97. Time required about 28 sec. 
987654321 = 3 * 3 * 17 * 17 * 379721. Time re- 


quired about 220 sec. 
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EUCLIDEAN ALGORITHM 


1. Purpose 


To determine the greatest common divisor of two inte- 


gers a and b, not both zero, |ja|, |b| < 10 < 


2. Method 


The Euclidean algorithm. We determine a sequence of 


nonnegative integers tn, } by 


Ny ?= max Clare bl). 

n, := min (jal,{bl]) , 

a ee _ a] 

eae _ = , 
k k-2 Ny kak 
k = 2, 3, ... . This sequence is decreasing, and thus 
it eventually reaches zero. If ny = 0, then 
e= nN 
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is the greatest common divisor of a and b (see ACCA II, 


§ 12.2). 


3. Flow Diagram 


max(|a|,|b|) 


:= min(|a|,|b]) 


Yes 


Display 


Beal 
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4, storage and Program 


Ry R, R, R, R, Re Ry 
[2] 
9 a 

00 25 GTO 00 
01 ABS 26 
02 xiy 27 
03 ABS 28 
04 x*“<y 29 
05 xiy 30 
06 sTo 0 31 
07 cae 32 
08 sTol as 

+09 RCL1 34 
10 x=0 35 
11 GTO 24 36 
12 RCL O 37 
13. RCL O 38 
14 RCL1 39 
15 + 40 
16 INT 41 
17. RCL 1 42 
18 sto 0 43 
19 x 44 
20 Z 45 
21 sto l 46 
22 PAUSE 47 
23 GTO 09 48 
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Ss Operating Instructions 


Load the program. Switch to RUN. Press 
FIX 0 


to get integers displayed as integers. Load a into 


Ro and b into Ri: Press 
PRGM 
R/S 


to start computation. The calculator will stop by 
displaying c, the greatest common divisor of a and Db. 


If both a and b are zero, c = 0. 


6. Examples and Timing 


a = 45, b = 96 yields c = 3. Computing time 
about 3 sec. ; 

[2 | a= - 965302379, b = 980051 yields c = 997, the 
correct result. Computing time about 5 sec. 

(3) Ae Lie b= O yields-c = 1, 
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RATIONAL BINOMIAL COEFFICIENTS 


l. Purpose 


To represent, for arbitrary rational pop = a the bino- 


mial coefficients 


_ /?P = 
by (e) = (0) Ly 
_ /f . PCO=1)) (0=2), 4.5 (0=ntl) 
n= 1], 2, ..., aS irreducible rational fractions, 
rn 
bo) Sy (1) 
n 
where r_ands_ are integers, (r_,sS_) = 1, and 
n n nn 


So > 0. Such rational values of the binomial coeffi- 
cients are often needed in analytical work involving 


the binomial series, 


(l+x)? = & DAG go Bl XS e 4 
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2. Method 
The recurrence relation By (0) = 1; 
_ pon — 
bia 8) = apy Pyle) » B= 0, 1, 2, eee 


is used. If 90 is a rational number, po = a and if be 


is represented in the form (1), this yields 


r Yr 
n+l _~p- qn on 
Snel q (n+l) s 


We thus first compute the integers 


* ee 2 x = 
cael ° (p qn)r. r S441? q(n+1)s_ ; 


then by means of the Euclidean algorithm determine 


their greatest common divisor 


* * 
ron (cael! eae 
and find 
* * 
es ee 5. = nel 
n+1 lc | ‘ n+l lc] 
The process is started with r, = s, = 1. The short 


0 0 
version of the Euclidean algorithm used here is bor- 


rowed from the program "Exact continued fractions for 
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quadratic irrationalities". 


3. Flow Diagram 


C= -(r*os*) 


(Euclidean algorithm) 
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4, Storage and Program 


Ro R, R, R3 Ra R, Re Ro 
r,r* s,s* n r* S 

00 25 STO 7 
O01 CLX > 26 RCL 6 
02 STO 4 27 RCL 6 
03 1 28 RCL 7 
04 STO 2 29 <7 
05 STO 3 30 INT 

> 06 RCL 2 31 RCL 7 
07 PAUSE 32 STO 6 
08 RCL 3 33 * 
09 R/S 34 7 
10 RCL 0 35 STO 7 
11 RCL 1 36 x # 0 
12 RCL 4 37 GTO 26 
13 = 38 RCL 6 
14 = 39 ABS 
15 STO*2 40 STO+2 
16 1 4l STO+3 
17 STO+4 42 GTO 06 
18 RCL 1 43 
19 RCL 4 44 
20 * 45 
21 STO*3 46 
22 RCL 2 47 
23 STO 6 48 
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5. Operating Instructions 


Load the program; move operating switch to RUN. Press 


to obtain integer display. Load numerator p of o into 


R, and denominator q into R On pressing 


0 1 
PRGM 
R/S 
calculator briefly displays ry = 1 and halts by dis- 
playing So = 1. Pressing 
R/S 


will cause display of r, and, after a short pause, 


1 


S etc. Pressing 


1’ 


after display of S. will cause the decimal value of 


K 


bn(p) = a 


yn) 


to be displayed, without disturbing subsequent compu- 
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tations. Also, after display of sat the current value 


of n may be exhibited by pressing 
RCL 4 


To re-inspect oe after display of Si, press 


m 
AV 
KK 


For large n, integer overflow may cause the values of 
ry and Sn to be inaccurate; however, the decimal va- 


lues of be) will still be accurate. 


6. Examples and timing 


(°) = Ly 6; 155-20, 155 65 Ly 05. 07 Us cas 
1 
p=5 , i.e p=1, q = 2 yields the values 
M2 Goo a. eh peat 2, ek 
n Der 8’ 16’ 128’ 256’ 1024’ 
33 _ 429 715 _ 2431 4199 
2048’ 32768’ 65536’ 262144’ 524288 
eee 52003 — 185725 
4194304’ 8388608’ 33554432’ 
334305 9694845 20036013 


67108864’ 2147483648’ 4867629602° 


Computing time approx. 165 sec. 
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CONTINUED FRACTION REPRESENTATIONS OF REAL NUMBERS 


dees Purpose 


Given a real number po, to find its standard represen- 


tation either by a terminating continued fraction, 


if p is rational, or by a nonterminating continued 


fraction, 


0 = b, + + PY tee. o> (2) 


if p is irrational. Here bo is an integer, and the 


b (i = 1, 2, ... ), if defined, are positive integers. 


2. Method 


By denoting by [x] the greatest integer not exceeding 
x, the be may be computed as follows (see ACCA II, 
§ 12.2): Let 
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by :_= [oe] v aa := p 7” Do a 
and for k = 1, 2, ... , if Py + 0, 
1 1 
b, := an Pan c= —-— - db, . (3) 
k Py k+1 Py k 


If p is an integer, then already Pp, = 0, and the re- 


presentation (1) reduces to 9 = b If pop is rational 


o° 
but not an integer, then eventually op = 0 for some 


n+1 
n> 0. The b b be thus generated are the 


0’ 1! eco fF 
partial denominators of the continued fraction (1). 


If 9 is irrational, then Py # 0 for all k 2 1. 

3. Flow Diagram 

A slight complication arises through the fact that 
the operation INTp agrees with [p] only if p 2 0, or 


if p < 0 and p is an integer. For nonintegral p < 0 


we have 
oy, = INE a 


Otherwise, the flow diagram is straightforward. 
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FRACo = 0 ? 


Display 
:= INTo 


Display 


bo := INTo - l 


Yes 


Display 
1 


Ne ee 
Pk 


Display 0 


4. 


Continued Fraction Representations 


Storage and Program 


R 


0 


Ry 


00 
Ol 
02 
03 
04 
05 
06 
07 
08 
09 
10 
11 
12 
13 
14 
15 
16 
17 
18 
19 
20 
21 
Z2 
23 
24 


27 
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Bs Operating Instructions 


Load the program; switch to RUN. Load op into the X 


register. Press 
FIX 0 
to get the integer display. By pressing 


PRGM 
R/S 


Do will be displayed. Pressing 


R/S 


repeatedly will display b b If a zero is 


; ge! Svinte 
displayed, this means a at algorithm has termina- 
ted, the number op being rational. Note, however, that 
because of rounding errors the algorithm does not al- 
ways operate as predicted by the theory. Thus for ra- 
tional p we may sometimes get an exceedingly large 


bol in place of 0, indicating that on the computer 


Pn+1 
irrational p, the b. from a certain point onward ge- 


was not exactly zero as it should have been. For 


nerally are wrong. These deficiencies could only be 
corrected by working with multiple precision (compli- 
cated on the HP 25) or, in certain special cases, 


with integer arithmetic (see the following program) . 
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6. Examples 

[1] p» = 2/3 yields by = 0, 1, 2, 0. Indeed, 
aie 
3 1 2 

p = -4.5 yields b, = -5, 2. Indeed, 
-4.5 = -5 +5 


p = 26/47 yields b, = 0, 1, 1, 4, 4, 1, 
(20000000). The last entry is wrong. In fact, 


26 1] i| 1 | 1 | 1 | 
“7 tro fa he eT 


[4 | For p = e we get, as discovered by Euler, ba = 


The last entry should be 10 and hence is wrong. 
o = 7 furnishes b, = 
Sy: “Ty boy hy. G29) 


The first few approximants of the continued frac- 


tion for 7 thus are 


30 
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a - 22 = 3.14285714 , 
7 7 
1 1 89S ox 

3: 4 + nH = 106 = 3.14150943 , 


1 di L 355 
3 +H ged 113 3.14159292 


The next approximant, unfortunately, is already 
wrong (b, should be 292), which is indicative of 
the great accuracy that is already achieved by 


the foregoing loworder approximants. 


 ~1+75 _ = _ a 
0 = 5 . We get bo = by =... boy = 1, and 
then b,5 = 2, which is wrong. 
0 = Y7. The sequence of bis is as follows: 


The last two examples illustrate the fact that 
the continued fractions representing quadratic 
irrationalities are ultimately periodic. The 

following program shows how to construct them 


without rounding error. 
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EXACT CONTINUED FRACTIONS FOR 
QUADRATIC IRRATIONALITIES 


Lg Purpose 


To construct, without rounding error, the standard 
continued fraction representation of a quadratic ir- 


rationality, that is, of a real number 


0 = p+ qvr (1) 


wherr p, q, mM, Yr are integers and r is not a square, 

m #0. The set of quadratic irrationalities coincides 
with the set of irrational solutions of quadratic 
equations ase + bx + c = 0 with integer coefficients 
a, b, c. By a celebrated theorem of Lagrange (see 

ACCA II, § 12.2) the standard continued fraction re- 
presentation of a quadratic irrationality is ultimate- 
ly periodic. It thus suffices to determine, in addi- 
tion to the nonperiodic part, one period of the con- 


tinued fraction. 
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2. Method 


The same division algorithm as in the preceding algo- 
rithm, but performed in rational arithmetic. If r is 
fixed, the numbers of the form (1) create a field F; 
that is, sums, differences, products, and quotients 
of such numbers can again be represented in the same 
way. This is clear for sums, differences, and pro- 


ducts; for quotients, we may use the identity 


_ m _ = mp + mqvr 
p + qvr q’r - p* 


Ole 


Each number p € F thus is represented by a triple of 
integers [P, qd, m]. Many triples represent the same 0, 
but among all triples representing op there is a unique 
representation, which we call the reduced representa- 
tion, that is characterized by the fact that p, q, m 
have no common factors and that m > 0. We may use the 


symbol 


[p, dq, mj] v 1 
for any representation of po, and 
[P, 4, m| % p 


for the reduced representation. For any extended nu- 


merical work in F one should use reduced representa- 
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tions, both to prevent overflow and to recognize re- 
peated elements. To this end, the integers p, q, m 
should be divided after each operation in F by their 
greatest common divisor (g.c.d.), say, d. This is de- 
termined by first calculating the g.c.d. (p,m) of p 
and m by the Euclidean algorithm, and then construc- 
ting d = ((p,m),q) by another application of the Eu- 
clidean algorithm. 


For our present purposes the division algorithm 


is formulated as follows: Let OL := 0, and for k = 0O, 
1, 2, .«... compute 
by := [o,] j 
O = - 
k+1 A by 


There is no possibility of vanishing denominators, be- 
cause the continued fraction is a priori known to be 
infinite. The period of the continued fraction is com- 


> : 
pleted as soon as Qin = Oy = 0 and n > O, that is, 
as a reduced triple representing Oy is repeated. 
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3. Flow Diagram 


Display 


b = [a] 


Inspect 


P, Gq, Mm 


Display 


(p,m) 


Display 


((p,m) ,q) 


To save programming space, the program computes by 
correctly only if po > 0. For similar reasons, a 
branching after the second application of the Eucli- 


dean algorithm must be performed manually. 
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4. Storage and Program 


Ro R, R, R Ra Re Re Ro 
p q m r temp temp temp 

00 25 R + 
Ol RCL 3 26 STO 2 
02 Yx 27 STO 7 
03 RCL 1 28 RCL 0 
04 ~ 29 STO 6 
05 RCL 0 + 30 RCL 6 
06 + 31 RCL 6 
07 RCL 2 32 RCL 7 
08 > 33 = 

09 INT 34 INT 
10 R/S 35 RCL 7 
11 RCL 2 36 STO 6 
12 = 37 = 

13 STO-0 38 = 

14 RCL 1 39 STO 7 
15 a 40 x # 0 
16 RCL 3 41 GTO 30 
17 * 42 RCL 6 
18 RCL 0 43 R/S 
19 er 44 RCL 1 
20 - 45 STO 7 
21 RCL 2 46 GTO 30 
22 STO*L 47 STO+0 
23 CHS 48 STOT1 
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A short version of the Euclidean algorithm is con- 


tained in instructions 30 through 43. 


oe Operating Instructions 


Load the program, then move the operating switch to 


RUN. Press 


FIX 0 


to get the integer display. Load 


p into Ro 
i R 
q into 1 
m into R, 
into R 
r into 3 


The following instructions should now be carried out 
cyclically. *Press 
PRGM 
R/S 


The computer very shortly will display b Pressing 


0° 
R/S 


will cause program to compute nonreduced representa- 


ti On) py a5, Jer 1/(a) - by) and to compute and dis- 
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play (p,m). Pressing 


R/S 


once more causes display of d = ((p,m),q). Now press 


GTO 47 
R/S 


This will compute the reduced representation [p,q,m|] 


% as: The calculator will stop by displaying d. To 


inspect reduced representation, press 


RCL 0 
RCL l 
RCL 2 


This will in turn show p, g, m. Now go back to * to 


compute b, and so forth. The period of the fraction 


1 
is complete if a reduced triple p,q,m is repeated; 
the period begins where the repeated triple occurs 


the first time. 


6. Examples and Timing 


0 = 24 — 715 . We obtain 
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b, 1 5 2 3 2 
1 

P, 24 7 3 | 3 

qs -1 1 1 1 

m 17 2 3 3 


The repeated triple is encased. Indicating the period 


by a bar, the sequence of the be thus is 


b, 0 i 8 6 7 1 1 
Ps 4 28 17 15 15 13 6 
q. i -7 7 7 7 7 7 
m, 7 11 4 5 r 19 11 
b, it 30 if 1 1 7 6 
P, 5 15 15 5 6 13 15 
q. 7 7 7 7 7 

20 1 20 i 19 4 5 


Total computing time about 10 min. 


Part 2 
ITERATION 


Tteration 41 


ITERATION 
4 Oe Purpose 


To determine a fixed point of a given function f, 


that is, a solution of the equation 
x = £(x). (1) 


This may also be used to construct solutions of equa- 


tions of the form 
g(x) = 0 (2) 
by seeking fixed points of the function 
f(x) := x + c g(x), 


where c is a suitably chosen constant or function. 


2. Method 


We construct an iteration sequence {x} by choosing 


Xo arbitrarily and forming 
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If f maps a closed interval I into itself, if f is 
contracting on I (i.e., if there exists k < 1 such 
that |£(x) - £(y) | < k|x - y| for arbitrary x, y € I), 
it can be shown that the equation (1) has precisely 
one solution s in I, and any iteration sequence star- 


ted with an Xp & I converges to s (see ENA, § 4.2). 


3. Flow Diagram 


n: n+1 
ot 3 f (x) 


Display x 


The index n merely counts the iteration steps; it is 


not actually needed in the computation. 


Iteration 43 


4, Storage and Program 


Ro Ry R, R, Ry Re Re Ro 
x nN 
00 08 PAUSE (R/S) 
Ol STO 0 09 GTO 04 
02 CLX 10 
03 STO 7 11 
04 1 12 
05 STO+7 13 
06 GTO 10 14 
07 STO 0 15 


The program to compute £ should be in the locations 
10 through 49. The last instruction of this program 
should be GTO 07. 


3 Operating Instructions 


Load the program; turn the switch to RUN. Select the 


mode of displaying numbers, for example, by pressing 


FIX 8 


Load the starting value x, into the X register. Pres- 


0 
Sing 
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PRGM 
R/S 


will start the computation. The computer will pause 
briefly after computation of each iterate and display 


it. If a stop is desired, instruction 08 should be 


R/S; in this case, 


R/S 


must be pressed after each display. The iteration 


index n is displayed by pressing 
RCL 7 
after display of x. No convergence test is provided 


by this simple program. 


6. Examples and Timing 


f(x) := /2 + x . The program to compute f is 
10 RCL 0 
11 2 
12 + 
13 Vx 


14 GTO 07 


Starting with Xo 


aac 
*14 
*15 


0, we get 


1.41421356 
1.84775907 
1.96157056 
1.99036945 
1.99759091 


1.99999996 
.99999999 
2.00000000 


Iteration 


which, to the number of digits shown, equals s 


= 2. Computing time about 20 sec. 


f(x) := I . ae Again starting with XQ = 0 we 
get 
x1 1.00000000 
0.50000000 
0 .66666667 
0 .60000000 
X19 0.61803400 
X59 0.61803398 
X54 0.61803399 
X59 0.61803399 


This equals s = 5 (v5 - 1), the unique positive 
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solution of 


Computing time about 30 sec. 


f(x) := (x - ae With x, := 3 we get a rapidly 
divergent sequence; with Xy t= 2 we find 
x = 1.0000 
Xy = 0.0000 
Xx. = 1.0000 
Xy = 0.0000 , 


an oscillating sequence. 


[4] £(x) := ee: X,) i= 0 produces a sequence that, 
although theoretically convergent, does not in 
fact converge. After n = 833 iterations, the 


sequence cycles between 


0.99954403 and 
= 0.00004561 


833 
*834 


Iteration with Aitken Acceleration 47 


ITERATION WITH AITKEN ACCELERATION 


1 Purpose 


To determine a fixed point of a given function f more 


rapidly than with ordinary iteration. 


2. Method 


Aitken acceleration. Along with the sequence tx}, 
which is generated as in the preceding program, we 


generate the sequence of accelerated values {xi}, 


where 
(Ax ie 
nes n 
KP r= kL 5 
A 
n 
(A = forward difference operator). If {x} converges, 


then under certain conditions (see ENA, § 4.4) {x} 


converges to the fixed point more rapidly than {x}. 
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3. Flow Diagram 


Display X, i= ¥ 


Compute and display Xo 


4. Storage and Program 


Iteration with Aitken Acceleration 


GTO 


49 
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Program locations 39 through 49 are reserved for the 
program computing f(x). This program should assume x 
in R, and should put f(x) into Re - Only the stack may 
be used for temporary storage. The last instruction 
should be GTO 06. The program for computing the f of 


example is shown above. 


5. Operating Instructions 


Load the program (including the program to compute f); 
Switch to RUN. Select the mode of displaying numbers, 


for instance 


FIX 8 


If £f contains trigonometric functions with argument 


in radians, press 


Load the starting value Xo into Ro - When 


PRGM 
R/S 


is pressed, the computation starts. The calculator 
will briefly display each X42 and stop at display 


of xO Press 
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R/S 


to start a new cycle. No convergence test is provided 


nor are iterations counted. 


6. Examples and Timing 


f(x) := /2 + x . Sequence {xi} .: ae Xy = 0: 
03942606 
-00208254 
00012537 
.00000776 
-00000048 
-00000003 
-00000000 


~ Mm Mm Mm mM Mm M 
A-M- B= We N-H- O- 
II 
Nm NO MO NH NHB HW WB 


The fixed point s = 2 has now been reached in 

Six iterations (which requires the computation 
of 8 iterates x): Computing time is about the 
same as before, because of the additional work 


required to compute x 


[2] £(x) := ee” . With x, = 0 we get 
XK. = 0.36787944 xi = 0.61269984 
0.69220063 0.58222610 


0.50047350 0.57170577 


52 Iteration 


At this 
whereas 


correct 


= 0.60624354 x4 = 0.56863881 
= 0.56706790 X14 = 0.56714330 
= 0.56718605 Xa = 0.56714329 


point the sequence {xi} has converged, 
the sequence {x} still has only four 


digits. Computing time about 80 sec. 
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AITKEN-STEFFENSEN ITERATION 


f- Purpose 


To compute the fixed points of smooth functions f 


without restrictions on the slope of f. 


2. Method 


Aitken-Steffensen iteration. This is a variant of the 


Aitken acceleration considered in the preceding pro- 


gram. Each x} is used as a starting point for a new 


0 
iteration. Equivalently, we generate a sequence 
ae by choosing oe arbitrarily and forming 
— (n), — _(np2 
(n+1) (n) EC eee wae 
= om fs ~ (n) (n) Gay = 
F(f(x °° )) - 2£(x ) + x 
(0) 


If f has a fixed point s, if f'(s) #1, and if x 
(n) 


is chosen sufficiently close to s, then x > Ss with 


quadratic convergence (see ENA, § 4.11). If f'(s) = 1, 


x) 


convergence is dubious; in particular, some may 


fail to be defined because of vanishing denominator. 


54 Iteration 


3. Flow Diagram 


Display x 


Aitken-Steffensen Iteration 


4, Storage and Program 


Xo 


Ps oh Rs Re 
AX, A*x, x 4 
25 RCL 3 
26 mG 
27 RCL 4 
28 = 
29 CHS 
30 RCL 0 
31 + 
32 STO 0 
33 R/S 


36 1 
37 0 
38 * 
39 CHS 
40 e* 
41 sTO 6 
42 GTO 06 
43 

44 

45 

46 

47 

48 
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The program to compute f£ should be in locations 35 


through 49. The value of x is in R y = £(x) should 


5? 
be put in Re - Use the stack only for temporary storage. 
The last instruction is GTO 06. The program for exam- 


ple [4] is shown above. 


5. Operating Instructions 


Similar to preceding program. After 


PRGM 
R/S 


is pressed, the computation starts. The computer dis- 


(n) 


plays each x and stops. Press 


R/S 


to start new cycle. No convergence test is provided. 


6. Examples and Timing 


fh: eee oe S00 yields sequence pee 4 


as follows: 


= 2.03942606 
= 2.00000802 
x3) = 2.00000000 


Aitken-Steffensen Iteration 5.7 


The limit has been reached (to within machine ac- 


curacy) in three iterations. Computing time 10 


Sec. 
1 (0) _ | 
f(x) := i x = 0 yields 
<'-? = @e6eeeees 
x2) = 0.61818182 


x3) = 0.61803399 


Computing time 15 sec. 


Lx) Sc x = i : oo = 3 yields 


(1) 


x = 2.75000000 
.63888889 
.61864329 
-61803453 


-61803399 


bh NO NO N LH 


The last value is = (3 + /5), the larger fixed 


point of f. Starting with oe 


= 0 we get 
.50000000 
.38888889 
. 38199234 
.38196601 
-38196601 


> a SP a > ae SE a 


The last value is =(3 - /5), the smaller fixed 


Iteration 


point. Computing time in both cases about 15 sec. 


£(x) := Soe With a0) = 0 we get 


x 4) = 0.50001135 
0.32882657 
0.23868709 
0.19120788 
0.17597093 
0.17456388 
0.17455280 
0.17455280 


Iteration converges after eight steps. Computing 


time about 30 sec. 
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NEWTON ITERATION FOR COMPLEX ROOTS 


1. Purpose 


To demonstrate the convergence of Newton's method to 
compute the square root of a given complex number 
c=a+t ib # 0. This method consists in forming the 


sequence {zo} = {x + iy} by choosing Zo and then 


computing recursively 


This sequence converges to the value of /c nearest to 


Zo if there is precisely one such value. It diverges 


if Zz) lies on the straight line through O perpendicu- 


lar to the straightline segment joining the two values 


of Yc (see ACCA I, § 6.12). 


2. Method 


We evaluate the sequence (1) until |z | < €, a 


n+l ~ 7n 
prescribed tolerance. The number of iteration steps is 


counted. 
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3. Flow Diagram 


Display Zo 
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ae: storage and Program 


a0 Ay RD ae Ra 3 Re a 
a b x y E | z ie Ax n 
0 0 n n 
00 25 RCL 2 
Ol CLX 26 * 
02 STO 7 27 RCL 0 
> 03 1 28 RCL 3 
04 STO+7 29 ~ 
05 RCL 3 30 = 
06 RCL 2 31 RCL 5 
07 + Pp 32 = 
08 0 33 RCL 3 
09 STO 5 34 - 
10 RCL 0 35 2 
11 RCL 2 36 = 
12 * 37 STO+3 
13 RCL l 38 RCL 6 
14 RCL 39 STO+2 
15 * 40 > p 
16 + 4l STO 6 
17 RCL 5 42 RCL 2 
18 = 43 PAUSE 
19 RCL 2 44 RCL 4 
20 - 45 RCL 6 
21 2 46 x = y 
22 = 47 GTO 03 
23 STO 6 48 RCL 3 
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5. Operating Instructions 


Load the program. Turn the operating switch to RUN. 


Select the mode of display numbers, for instance 
SCI 8 

Load the data as follows: 
a into 


into 


x into 


YA DD DW 
Ww NY FY © 


Yo into 


Load € into Rae For instance, for € = co press 


When 


PRGM 
R/S 


is pressed, the calculator starts iterating, pausing 
briefly after each cycle while displaying 7 a After 


convergence has been achieved, the calculator will 


Newton Iteration for Complex Roots 


stop and display final xX To find Y,7 press 


AN 


To display number of iterations required, press 

RCL 7 
If all iterates are to be recorded, instruction 43 
should be changed to R/S. The calculator will then 
stop after each cycle, displaying x To display Ye 
press 


RCL 3 


To compute next iterate, press 


R/S 
6. Examples and Timing 
ee c= 3 + 4i, Zo =l, €= 16-7. The successive 
iterates are 

2.00000000 + 2.00000000 i 
1.87500000 + 1.12500000 i 
1.99632353 + 0.99387255 i 
1.99999968 + 1.00001144 i 
2.00000000 + 1.00000000 i 
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Iteration 


Convergence has been achieved. Computing time 


35 sec. 
c and € as before, Zz = -i. Algorithm converges 
to 
- 2.00000000 - 1.00000000 i 
in n = 4 iterations. Computing time 40 sec. 
c = -1 + 0.00001 i, 2. = 1, € = 10°’. Algorithm 


0 
converges to 


0.00000500 + 1.00000000 i 


in n = 24 iterations. Computing time about 2 min. 


c = -l, z, = 1. No convergence. 


Part 3 
POLYNOMIALS 
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HORNER ALGORITHM 


i Be Purpose 


Given an arbitrary polynomial of degree S 4 with real 


coefficients 


p(x) = a,x + a,x +a 


and given an arbitrary real number x to determine 


0’ 
the coefficients b in the representation of p in 


powers of h := x - Xoo 


4 3 2 
P(X + h) = boh + bh + boh + b,h + Dy ; 


that is, the Taylor coefficients of p at Xo 


2. Method 


The Horner algorithm (see ENA § 3.4; ACCA I § 6.1). 
(m) 


n 


We generate coefficients b as follows. Let 
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and for m= 0, l, ... , 4: 


iM) ._ pim-1), 


(m) (m) (m-1) = 7 
0 e 0 g bo .— Xb 1 + bo w nN = DLs egg a m 
Then 
4-m 
bi = be eS Oe ik Sea. ap ad 


3. Flow Diagram 


Because an address modification is not available, the 
programming here is different from what it would be 
on a larger computer. By a cyclic permutation (here 
called rotation) of the avs the coefficients operated 
on are always found in the same registers. After the 
algorithm is executed, the be are stored where, pre- 


viously, the a, were stored. 
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Yes 


Display bo 


Stop 
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4. Storage and Program 


Ro Ry R, R Ry Re Re Ro 
ay ay a. a. ay Xo temp n.m 
00 25 RCL 0 
O1 CLX 26 STO 6 
02 STO 7 27 RCL l 
+ 03 RCL 7 28 STO 0 
04 INT 29 RCL 2 
05 STO 7 30 STO 1 
+ 06 RCL 7 31 RCL 3 
07 FRAC 32 STO 2 
08 9 3.3 RCL 4 
09 = 34 STO 3 
10 RCL 7 35 RCL 6 
11 + 36 STO 4 

12 3 37 

13 x < y 38 1 

14 GTO 25 39 STO +7 
15 RCL 0 40 

16 RCL 5 41 5 

17 * 42 RCL 7 

18 RCL 1 43 FRAC 

19 + 44 x <y 

20 STO 1 45 GTO 06 
21 3 46 1 

22 RCL 7 47 STO +7 
23 x2 48 GTO 03 
24 GTO 49 > 49 RCL 0 


Note: A fractional index n.m is used to save storage. 
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5. Operating Instructions 


Load the program and switch to RUN. Load the coeffi- 


cients an into Ro (n = 0, l, ..., 4) and Xo into Res 
When 
PRGM 
R/S 


is pressed, the calculator computes the ba and stops 


by displaying b At the same time, the DT are stored 


0° 
in Ran (m = 0, 1, ..., 4), thus enabling the operator 
to continue the computation immediately with a dif- 


ferent Xo The original a are lost. 


6. Examples and Timing 


Consider the polynomial 


p(x) = =. - A? + 3x- - 2x + 5 


Carrying through the algorithm with Xy = 0.5 


yields the Taylor coefficients at 0.5, 
b = 1.0000, -2.0000, -1.5000, -1.5000, 4.3125 
(time required about 29 sec). The last coeffi- 


cient is p(0.5). Repeating the algorithm with 
the foregoing data yields the Taylor coefficients 
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at l, 
bd = 1.0000, 0.0000, -3.0000, -4.0000, 3.0000 
Repeating three times with x, = -0.33333333 


0 
yields the coefficient arrays 


1.0000, SL e3 353% SLi SIoo 7 =2,.1461, 4.0123 
1.0000, -2.6667, =04. 393935) -1.1852, 4.5309 
1.0000, -4.0000, 3.0000, -2.0000, 5.0000 


The last array is identical with the initial 
array of coefficients, demonstrating the "group 


property" of the Horner algorithm. 


For the polynomial p(x) = x*, XQ = 1 we get 


b = 1.0000, 4.0000, 6.0000, 4.0000, 1.0000, 


the binomial coefficients (i Computing time 


29 sec. 
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NEWTON'S METHOD FOR POLYNOMIALS 


1. Purpose 


To determine all real zeros of a real polynomial of 


degree 4, 


p(x) = a,x +a +a e + a,x ta,. (1) 


2. Method 


Newton's method, using Horner's scheme to evaluate p 
and p' and deflating the polynomial after each zero 
has been found. Newton's method requires forming the 


sequence {x} according to 


p(x) 
x =x oo (2) 
n+l n p'(x_) ’ 
n 
with an arbitrarily chosen starting value x The se- 


0° 
quence {x} will converge to a zero, provided a real 
zero exists and XQ is chosen close enough to the zero. 


[Unless p'(0) = 0, the choice x, = 0 will frequently 


0 
work and will produce convergence to the zero of 
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smallest modulus. | For the purpose of this algorithm, 


Horner's scheme may be formulated as follows: Compute 


the array tb. } by 


b_, ?= 0, b, i= xb + dee go OK Sy. A 24°35 43 


then compute the array tc. } by 


Then p(x) = Dy and p'(x) = C.- The iteration (2) is 


terminated if 


_ -6 
Ib,| = |p(x)| < 10°. 


(This tolerance may be changed to any one-digit power 
of 10, but too stiff a tolerance may cause the se- 
quence {x} to cycle instead of to converge.) If 


p(z) = 0, the be are simply the coefficients of the 
deflated polynomial 


.. P(X) _ 3 2 
P, (x) oe box + b x + box + b, 


Because the coefficients b. are available, the de- 
flation thus is achieved by moving the b. into the 


locations of the a. (see operating instructions). 


Newton's Method for Polynomials 75 


3. Flow Diagram 


Compute p(x) 


Yes 


Compute p' (x) 


Display 


Display x 
Stop 


The program presents a particular challenge because 
there are ten essential quantities a 


b b b b 


or Syn Soe Aas 


Agr 1’ Por Par Bye x that have to be saved. 
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4. Storage and Program 


IV OV 


48 
49 
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Be Operating Instructions 


Load the program. Turn the operating switch to RUN. 
Select the mode of displaying numbers, for instance, 


by pressing 
FIX 8 


Load the coefficients a. of given polynomial into 
Rog = O, ... , 4. [Be sure to use indexing of co- 


efficients as in (1).] Load the starting value Xo 


into the X register. Press 


PRGM 
R/S 


to start computation. The iterates x will be compu- 
ted; each x will be displayed briefly. The calcula- 
tor will stop if x, meets convergence test |p (x) | < 
io: displaying final value of xO: Exponent -6 in 
convergence test may be changed to -m by changing in- 
struction 23 to m, where m is any one-digit integer. 
To deflate polynomial when zero has been found, 


press 


RCL 0 
STO 1 
RCL 5 
STO 2 
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RCL 6 
STO 3 
RCL 7 
STO 4 
CLX 
STO 0 


Load the starting value x, into the X register (if 


0 


x, = 0, it is already there) and press 


0 
R/S 
to restart computation. The process may be repeated 


until all real zeros have been found. 


6. Examples and Timing 


p(x) := = = ies: + 9557 - 96x + 24. 
Starting with Xy = O we get the smallest zero 
Zz, = 0.32254769 


after four iterations. Deflating and again star- 


ting with Xo = 0 yields 


Zz. = 1.74576110 


after five iterations. Deflating and starting 
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with XQ = O after five iterations yields 


Z. = 4.53662030 


Deflating once more yields 


9.39507092 


N 
II 


in one iteration. Computing time about 4 sec 


per iteration. 


p(x) := - = 10x> + 5° - 50x + 24 
= (x - 4) (x - 3) (x - 2) (x - 1). 


The program correctly finds the zeros 


Zz, = 1.00000000 
Zo = 2.00000000 
Zz. = 3.00000000 
ea = 4.00000000 


Total elapsed time (including deflations) 


2.75 min. 
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BERNOULLI'S METHOD FOR SINGLE DOMINANT ZERO 


l. Purpose 


To determine the zero of greatest absolute value of a 


real polynomial of degree 4, 
4 2 
p(x) = X + a,x + a,x + a.xta, , (1) 


if there is a single such zero. 


2. Method 


Bernoulli's method (see ENA, § 7.1; ACCA I, § 7.4). 


A sequence {x} is generated by choosing starting 


values Xyr Xp Xr KX, "arbitrarily" and letting 
x = = (ajx_y + aX 9 + a,X_3 + A,X _4) : 
n= 4, 5, ... . One forms the guotients 
xn 
DA. ; (2) 
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if p has a single dominant zero w and if the starting 
values are not in a certain set of measure zero, then 
the q, are defined, at least for sufficiently large n, 


and 


This may cause accidental termination if two consecu- 
tive q's happen to be identical. The program will al- 
so halt if an x becomes accidentally zero, or if 
there is overflow or underflow because x exceeds the 
range of the computer. 

Experience has shown that for polynomials with 
small integral coefficients these accidents most often 
are caused by a too special choice of the starting va- 
lues, such as Xy = X) = Xo = X3 = 1. The program, 
therefore, generates its own highly irrational star- 


ting values. 


3. Flow Diagram 


A difficulty is generated by the fact that all avai- 
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lable storage locations are taken up by the four co- 
efficients of the polynomial and by the four currently 
needed x: Therefore, there is no orderly transition 
from n to n+ 1; rather, each x is overwritten with 


x as soon as it is no longer needed. Furthermore, 


n+l 


the quotients q.. and qna-1) @re recalculated each time 


they are needed. 


Generate x x 


OF als 


Display w=q 


2 


4. 


storage and Program 


+ 09 
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45 
46 
47 
48 
49 


CHS 
STO 4 
RCL 5 
PAUSE (or R/S) 
ENTER 
ENTER 
RCL 5 
RCL 6 


GTO 00 


84 Polynomials 


5. Operating Instructions 


Load the program, switch to RUN, and choose mode of 


displaying numbers desired, for instance, by pressing 
SCI 8 


to get floating eight-digit display. Another possibi- 
lity is initially to display few digits only, which 
makes it easier to check convergence "by eye." Be sure 
to define coefficients of polynomial as in (1). Load 
the coefficient ay into Ry (k = 1, 2, 3, 4). Press 
PRGM 
R/S 


to start computation. The calculator will display 
briefly each qn: (If an indeterminate stop is desired, 
instruction 31 should be changed to R/S. In this case, 
R/S must be pressed after each display to continue 
computation.) The calculation will terminate if con- 
dition (3) is satisfied; the last qn will be dis- 
played. 

If starting values other than those provided by 
the program are desired, these should be loaded as 
follows: 

x into R 
into R 
Xo into R 


X3 into Ry 
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In this case, the calculation is started by pressing 


GTO 09 
R/S 
6. Examples and Timing 
p(x) := ra - LOK + so° - 50x + 24 


(x = 4)(x — 3) (x = 2) (xe = 1). 
The calculator displays 


w = 4.0000000 


after 60 iterations. Computing time approximately 
4 min. 

pe + a = ee - 15x + 18 

(x + ay a - 1)(x - 2). 


p(x) 


Convergence is very slow because the dominant 
zero has multiplicity > 1. Faster convergence is 
achieved (see ENA, § 7.4) by choosing the star- 


ting values 


Xy = 7 ay = = .3 

x, = 7 (2a, a5 a, Xo) = 23 

Xo = - (3a, + a5Xo + a,%,) = -—- 45 

X, = (4a, + A,X + ayX) + a,X>) = 179 


This yields 
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WS = 249999999 


after 45 iterations. Computing time approximately 


3 min. 
p(x) := rs = oe = on - 2x - 3 yields 


w = 3.6652758 


Computing time 3.4 min. 
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BERNOULLI'S METHOD FOR COMPLEX 
CONJUGATE DOMINANT ZEROS 


Le, Purpose 


To determine the quadratic factor formed by the two 
zeros of largest modulus of a real polynomial of de- 


gree 4, 


p(x) = a + a,x t+a ae + a.x ta, , (1) 


provided that all remaining zeros have smaller modulus. 


2. Method 


Bernoulli's method (see ENA, § 7.5), which in this 

Situation is identical with a pristine version of the 
quotient-difference algorithm (see ACCA I, § 7.6). We 
generate the sequence {x} as solution of the diffe- 


rence equation 
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n = 4, using arbitrarily chosen initial values Xr Xo 
Xor Xa. With the x we form the quotients 
x 
q n, n 
a ’ 
n X-1 


as in the preceding algorithm, and with the quotients 


the differences 


ae = Ty-1 "7 Gn-1 n aa 4n-29n-1 


If p has two dominant zeros w, and w and if all 


1 2" 
other zeros have smaller modulus, then the limits 


r := lim oe ; s := lims 
n-o n-©o 


exist and satisfy 
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in other words, Wy and Ww, may be determined as solu- 


tions of the quadratic equation 


Ox 3 (2) 


2 
x = xX aes 


Our program generates the values ra and Sr but 
it does not check for convergence because of lack of 
programming space. For the same reason, the starting 


x xX. muSt be supplied by the user. 


values x Xo 3 


Qo rat 


3. Flow Diagram 


Again, there is the difficulty that all available 


storage is taken up by the coefficients a, and by the 


k 
currently needed Xo: Auxiliary quantities, if they 

cannot be stored in the stack, are recalculated each 
time they are needed. The x, are overwritten with X41 


as soon as they are no longer required. 
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Display r := qs + qo 


Display s := 459) 
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4. Storage and Program 


00 25 = 

> OL. RCL 3 26 = 
02 RCL 7 27 RCL 5 
03 * 28 RCL 6 
04 RCL 2 29 4 
05 RCL 6 30 RCL 6 
06 STO 7 31 RCL 7 
07 x 32 = 
08 + 33 = 
09 RCL 1 34 = 
10 RCL 5 35 RCL 5 
iuae STO 6 36 RCL 6 
12 * 37 = 
13 + 38 x 
14 RCL 0 39 LAST X 
15 RCL 4 40 x < y 
16 STO 5 41 + 
17 % 42 PAUSE 
18 ae 43 LAST X 
19 CHS 44 RCL 6 
20 STO 4 45 RCL 7 
21 RCL 5 46 = 
22 = 47 2 
23 RCL 5 48 R/S 
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5 Operating Instructions 


Load the program; switch to RUN. Select the mode of 


displaying numbers (e.g., SCI 8). Load coefficients 
ay of given polynomial into R._4 (k = 1, ... , 4), 
being sure to use notation (1). Load irrational star- 
ting values (such as /t, e', V2, etc.) into Ryrosee se 
R.- Press 
PRGM 
R/S 


to start the computation. The calculator will alter- 
nate between a brief display of Ee and an indetermi- 


nate display of So: After each display of Sr 
R/S 


must be pressed to start new round of computation. 
Convergence of the sequences ix} and {s3 must be 
tested by eye. 

If indeterminate displays of both ry and Ss are 
desired, instruction 42 must be changed to R/S. If 
brief displays of both ry and S, are desired, instruc- 
tion 48 must be changed to PAUSE. 

An error halt may occur if an xX, or ane, becomes 
accidentally zero, or if there is overflow in x? In 
the first case, the calculation should be repeated, 
using different starting values. In case of an over- 


flow, the exponents in all x. should be scaled down 
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by the same number. 


6. 


Examples and Timing 


Let p(x) = er + 3x + LOK" - Vise“ + 7. We first 


use the display mode FIX 2 (for easier reading) 
and set both instructions 42 and 48 to R/S. Using 


the starting values 


XK, 3= - Sa Ry 
xX, i= Oe > Re 
= TT 
x =e >+ R 
1 Ys 6 
xX :=tane’" + R 
0 7 
we get 
xr Ss 
n n 
15.18 4.04 
= F405 = 354-37 
- 3.43 1132 
- 3.85 de eck: 
- 4.00 L3223 


Because convergence appears to be slow, we change 
both instructions 42 and 48 to PAUSE and let the 
computation go on. After a while, the sequences 
ir} and {so} will have converged to within 10." 
We change display to FIX 5 and continue. After 


convergence to within io? has been achieved, we 
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change display to FIX 8. After 
Es? 3.97856125 , Ss, = 13.36969365 


has been reached, values do not change anymore. 
Dominant zeros of polynomial thus are solutions 
of 


ee + 3.97856125 x + 13.36969365 


| 
ro) 
= 


that is, 


- 1.98928063 + 3.06797266 i 
Wigs SE 1.98928063 - 3.06797266 i 


_ _4 4 3 a2 2 1 
p(x) := x + EX + EX tex te 


After many iterations and about 18 min computing 


time, we get 


r = 0.2756646 , s = 0.4788911 
n n 


Thus p has the approximate quadratic factor 


7 - 0.2756646 x + 0.4788911 


with the zeros 


W, 5 = 0.1378323 © 66781544074 
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QUOTIENT-DIFFERENCE ALGORITHM 


1. Purpose 
To compute the zeros of a real polynomial of degree 4, 


p(x) = a,x + a,x +a i et a OR as (1) 


2. Method 


The progressive version of the qd algorithm (see ENA, 
§ 8.5; ACCA I, § 7.6). For the present purpose this 
may be described as follows: One computes a sequence 


of arrays of numbers 


by the initial conditions 


a 
(0)  . Jl (O) _ = = ; 
qy = ap ’ q5 = 92 = qd, = 0 ; (2) 
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and by the continuation rules 


nial = + a 4 Ee /k=1, 2, 3, 4, (4) 
where always a = ao := 0, and 
ge 
a := Ne ay eS. oye (5) 
dy 


(Note: These formulas are most easily remembered as 


"rhombus rules": 


The system of indexing used above is not identical 
with the indexing used in the sources cited.) 

Let the zeros of p be denoted by we und numbered 
such that 


IV 


> 
Jw | = [wo 


Then under certain conditions - for instance if the 
w, are all positive and distinct - all arrays a, 


exist, and 


Quotient-Difference Algorithm 97 
(n) 
lim q. = Wy ; 
n> 
k = 1, 2, 3, 4. Convergence of ie to Wh takes place 
n z 
with an error O(e€ ), where € := max(|w,.4/w |, 


Jw, /wy_4|)- By means 
also possible to use 


of complex conjugate 


3. Flow Diagram 


qs and e. by (2) and (3) 


Compute new qe by (4) 


Compute new e. by (5) 


Display n :=n+1 


of auxiliary computations it is 


the arrays an to compute pairs 


zeros (see ACCA I, § 7.9). 


Initialize 


98 


4. 


Polynomials 


Storage and Program 


0 Sa ae 
0 7 “2 
Ty sa 
00 
01 RCL 4 
02 RCL 3 
03 2 
04 STO 6 
05 RCL 3 
06 RCL 2 
07 * 
08 sto 4 
09 RCL1 
10 sTot2 
11 RCL O 
12 CHS 
13. sTotl 
14 CLX 
15 sTo 0 
16 STO 3 
17. sTo 5 
18 sto 7 
+19 RCL 2 
20 sTotrl 
21 RCL 4 
22 RCL 2 
23 - 
24 sT0+3 
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5. Operating Instructions 


Load the program, and move the operating switch to 
RUN. Select the mode of displaying numbers, for in- 
stance by pressing 


FIX 8 


(adequate here because of the convergence of the q's). 


Load coefficient ay into Ry Ck: = Op iby san -¢ 4). Be 
sure to define coefficients as in (1). Pressing 
R/S 


will start computation. The calculator will pause 
briefly after each array is computed, displaying its 


index n. To inspect array, press 

R/S 
during pause. The elements of the array computed last 
are then shown by pressing RCL k (k = ll, ..., 7). 
Pressing 


R/S 


will continue the computation. 
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6. Examples and Timing 


p(x) = x? - 16x? + 72x” - 96x + 24, the Laguerre 
polynomial of order 4. The initial array (set up 


by the program) is 


(n), 


The program produces the following values q. 


ae a” a ee 


10 9 .39729639 4.53453197 1.74562397 0.32254767 
20 9 .39507245 4.53661877 1.74576109 0.32254769 
60 9 .39507092 4.53662030 1.74576110 0.32254769 


Time per array (including pause) approx. 2.75 sec. 
A 3 2 
p(x) := x - 10x + 35x - 50x + 24 


= xe A CK ey TD) 


Values produced: 


- — a a a 

10 4.0759228 2.9543419 1.9711691 0.9985662 
20 4.0040290 2.9965044 1.99946 80 0.9999986 
30 4.0002261 2.9997832 1.9999908 1.0000000 
40 4.0000127 2.9999874 149999999 1.0000000 
50 4.0000007 2.999999 3 2.0000000 1.0000000 
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ROUTH ALGORITHM 


1. Purpose 
To determine whether a given real polynomial 
2 

Dix): 2= as f a. 4 aak” HP. eas ax 
of degree n s 7 is stable, that is, whether all its 
zeros have negative real parts. 
2. Method 
The Routh algorithm (ACCA I, § 6.7; ACCA II, § 12.7). 


For the purpose of this program, this algorithm may be 


described as follows. Let 


(0) 
by := ao, ; 
CI)... 
oj mes 
Ke Ou thy Ses. “CIE TS a4 Se = aa := 0.) For 


Sly 2y ate 4-3 do 
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4) 
Se ees, ee ok - 
4° asap Pk = 4G ka Pict! 
0 


The polynomial p is stable if and only if the numbers 


qe qo oe ey qn all exist and are positive. 


3. Flow Diagram 


i ec ea) 
k ; by = by ; 
storage space, quantities by are overwritten with by 


as soon as they are no longer needed. 


We set gq := Gye by := b To save 
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4. Storage and Program 


10 ms 35 
11 RCL 3 36 
12 STO 2 37 
13 = 38 
14 STO 1 39 
15 R ¥+ 40 
16 RCL 4 4l 
17 s 42 
18 RCL 5 43 
19 STO 4 44 
20 = 45 
21 STO 3 46 
22 R ¥ 47 
23 RCL 6 48 
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5. Operating Instructions 


After loading the program, switch to RUN. Load coeffi- 
cients a. of given polynomial into the storage regi- 
ster R, (i = 0, l, ... , 7). If degree n is less than 
7, the remaining storage registers are to be filled 


with O. Press 


PRGM 
R/S 


The computer will calculate and stop by displaying QW: 


Pressing 
R/S 


again will result in displaying Gos and so forth. The 
polynomial is stable only if the first n displayed qe 
are positive. The polynomial is unstable if some qe 

< 0 for i = n, or if display "Error" indicates divi- 


Sion by zero. 


6. Examples and Timing 


p(x) = 15 + 22x + ee + oe + ae The first 


four q, are as follows: 
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Os. 9 
0.08 


They are positive; the polynomial thus is stable. 


Indeed, its zeros are z = -l1 + iY2 and z = -2 +i. 


p(x) =1l+x + ie + oe + — + e + 2 + x . The 


following output is generated: 


1.00 
0.00 


Error 


Already q, = 0, indicating that p is unstable. 
p(x) =6 +x + re + rs + ae. + a + a + re 


The following q,'s are obtained: 


0.17 

= 10% 33 
- 24.50 
25.87 

= 0437 
= oO ele 
0.11 
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Not all q, are positive, indicating that p is un- 
stable. Because all q. exist, the fact that four 
q, are negative permits the conclusion that p has 
precisely four zeros with positive real part 
(ACCA I, § 6.7). 

Total execution time for each of these exam- 


ples is less than 20 sec. 
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SCHUR-COHN ALGORITHM I 


1. Purpose 
To determine whether or not all zeros Za Zor Zar Za 
of a given polynomial of degree 4, 
4 3 2 
= + 
p(x) aX + a,x + ax + aX ay, (1) 


satisfy |Z. | Se Ts 


2. Method 


The Schur-Cohn algorithm (see ACCA I, § 6.8), in 
slightly modified form for purposes of efficient com- 


putation. For a polynomial of degree n we would form 


a triangular array of numbers a Ci Og: dg et waie 2 TS 
k = 0, l, ... , n-m) as follows: Let 
BO sed WG By ane od 


Then for m= 0, l, ... , n-l do 


108 Polynomials 


, KO HO, Ly aexk 9) Me=meL. 


All zeros of p satisfy |Z, | < 1 if and only if the q_ 


exist and satisfy 


3. Flow Diagram 


The program implements the algorithm by means of a 


Onedimensional array of quantities b, (k = 0, l, ... , 


k 
n) which in case of a polynomial of degree n would be 


defined as follows: 


In terms of the b, the algorithm performs as follows: 


k 
For m= 0, l, ... , n-l do’ 
an 
Se ge ge = Aa Ae i 
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b/ := by - qb gp R= Te 2s , M3 (3a) 
by := by - cL aey er , K=am1, 6° As (3b) 
by = bred oe = he ee Pe (4) 


To carry out (3), the last n-l elements of the se- 
quence bi bos bee, oy bel! bos ele bor where be 
occurs m-l times, are stored in the stack of the cal- 
culator. On the HP-25, this is possible for n = 4. The 


program is redundant inasmuch by is calculated m+l 


times in the m-th cycle. 
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Put bo into stack 


Calculate b' by 


k 


No 


yes 


Show n, Stop 
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4. Storage and Program 


Ry Ry Ry R, R, R, Re Ra 
“6 “1 aD sae a4 

00 25. RCL 5 
01 CLX 26 * 
02 sto 7 27. +sTO-1 

+ 03 STO 6 28 R + 
04 RCL 4 29 RCL 5 
05 RCL O 30 x 
06 2 31 sTto-2 
07 sto 5 32 Ry 
08 PAUSE 33. RCL 5 
09 RCL 4 34 * 
10 * 35 sTo-3 
11  sfo-0 36 =o RCL 3 
ia - ‘Ren a 37. «STO 4 
is “ReE-o 38 RCL 2 
14 RCL 3 39 «=6sToO 3 

+15 RCL 6 40 RCL 1 
16 x= 0 41 sto 2 
17. GTO 24 42 RCL 0 
18 R 4 43 sto 1 
19 RCL 4 44 1 
20 1 45 STO+7 
21 sTo-6 46 4 
22 R ¥ AG Rem. 7 
23 GTO 15 48 x<y 


24 R ¥ AQ GTO 03 
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5. Operating Instructions 


Load the program. Move the operating switch to RUN. 
Select the mode of displaying numbers. (Because it 


only matters whether lq,] <1 or 1a | 2 1, FIX 2 is 


adequate except in doubtful cases.) Load coefficient 
ay of given polynomial (1) into Ry (Kose. Dy aren... g Ads 
When 
PRGM 
R/S 


is pressed, the calculator will carry out the algo- 
rithm, briefly displaying each qd. (k = 1, 2, 3, 4), 
and stop by showing n = 4. The zeros of p all lie in- 
side the unit circle if and only if all qx exist (no 
error halt due to division by 0), and if they satisfy 
la, | < 1. 

If it is desired to determine the number of zeros 
satisfying | 2,, | > 1, the quantities q, should be re- 
corded and used as input for the program "Schur-Cohn 
Algorithm II". Instruction 08 should be changed to R/S 
in this case, and R/S should be pressed after display 


of each q.- 


6. Examples and Timing 
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The algorithm yields, using FIX 1 display for 


quick reading, 


Computing time 14 sec. All a, | < 1, thus all 
zeros of p satisfy |Z, | < 1. This could also have 


been inferred from the Gauss¢lucas theorem (ACCA 
Ge ee 

I, § 6.5), since p(x) = ee a or from the 

Enestr6m-Kakeya theorem (ACCA I, § 6.4), because 

the coefficients a, are positive and from a mono- 


tonically decreasing sequence. 


Sea Ge oe a ae 


Using FIX 1 display the program obtains 


The value 1.0 is doubtful. We repeat the calcu- 
lation with FIX 8 display, obtaining 


0.25000000 0.33333333 0.50000000 1.00000000 


Not all lay, | < 1, hence not all zeros satisfy 


|Z, | < 1. Indeed, p(-1l) = 0. 
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SCHUR-COHN ALGORITHM II 


i Ee Purpose 


To determine the number of zeros Zs of a polynomial of 


degree 4, 


p(x) = ax + a,x +4 2 +a.xta, , (1) 


that satisfy |Z, | > 1, that is, that lie outside the 


unit circle. 


2. Method 


The Schur-Cohn algorithm (see ACCA I, § 6.8). Fora 
polynomial of degree n this runs as follows: If all 
numbers qs determined by the Schur-Cohn algorithm I 
exist and satisfy 1a, | #1, the number m of zeros out- 
Side the unit circle equals mos where m is determined 


recursively as follows: 
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Megat Oe eerie eo 
ete oe WA eee 


If some ay, | = 1, the number m cannot be determined by 


this algorithm. Our program then displays m = 99. 


3. Flow Diagram 


This is a straightforward realization of the formulas 

(2). Because the q, under examination always has to be 
at the same place, the remaining q, are shifted up one 
place at the end of each cycle. The program implements 


the case n = 4. 
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Display "99" yes 


Stop 


k = n? 
Yes 
Display m 


Stop 
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4. Storage and Program 


Ro Ry R, R, R, Re Re Ra 
4) To 43 Tq ‘ - 
00 25 sto 2 
01 CLX 26 GTO 04 
02 STO 6 > 27 9 
03 sto 7 28 9 

> 04 1 29 GTO 00 
05 STO+6 > 30 RCL 7 
06 RCL 4 31 GTO 00 
07 ABS 32 
08 x < Yy 33 
09 GTO 16 34 
10 x= 35 
11 GTO 27 36 
12. RCL 6 37 
13. RCL 38 
14 2 39 
15 sto 7 40 

+ 16 RCL 6 41 
17 4 42 
18 x = y 43 
19 GTO 30 44 
20 RCL 3 45 
21 sto 4 46 
22 RCL 2 47 
23. ~STO 3 48 
24 RCL1 49 
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pee Operating Instructions 


Load the program and move operating switch to RUN. 


Press 
FIX 0O 


for integer display of number of zeros. Load qe into 
Ris i=1, 2, 3, 4. When 
PRGM 
R/S 


is pressed, the calculator after a short time will 
display the number m of zeros outside the unit circle. 
If some 1a, | = 1 and hence the algorithm fails, the 


integer 99 will be displayed. 


6. Examples and Timing 


With the data of example of the preceding 


algorithm we get 


Indeed, we already know that this polynomial has 
no zeros outside the unit circle. Computing time 


about 4 sec. 
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p(x) = oe = 18x> + 52x + 3x. 15 
The program "Schur-Cohn Algorithm I" yields the 


values 
aq, = =0.6, -0.4, 4.7, 0.2 
Using these data, the present algorithm yields 


m= 2, indicating that precisely two zeros of p 


satisfy |Z, | Shs 


Part 4 
POWER SERIES 
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RECIPROCAL POWER SERIES 


1. Purpose 
Given a formal power series with leading coefficient l, 


P=l1+t+asxt+a a + 
—— a 2 eee v 
to compute the coefficients b of the reciprocal power 


series 


satisfying pee PS! le 


2. Method 


This is one instance where the very limitations of the 
pocket calculator make it necessary to invent a new 
and unorthodox algorithm to solve the problem. Ordina- 
rily one would use the fact that 


2 


(l + a,x ta ee + ...)(1 + Dlx + b,x ee “se ise e sees aly 


1 2 1 
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which implies for n= 1, 2, 
bo + aba + ab. _» + + ois 0:4 
to obtain the recurrence relation by = 1, 
Pag eae > “So aso = ~ ayPo » 
n= 1, 2, ... . However, to compute b by this method 


requires computing and storing all b, where k <n, 


which exceeds the storage capacity e the pocket cal- 
culator except for trivial values of n. The same is 
true for "fast" algorithms for computing the recipro- 
cal series uSing Newton's method and the Fast Fourier 
Transform (see Aho, Hopcroft, and Ullmann, The design 
and analysis of computer algorithms, Addison-Wesley, 
1974). 

Instead, we use the following approach which, how- 


ever slow the resulting algorithm, works on the HP-25. 


Writing 
Q:= a,x + a,x +a an + ; 
1 2 3 
so that P = 1 + Q, we obviously have 
pr ii Q + 9° = 9? Es 


For any n 1, the coefficient bo thus is the sum of 
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all coefficients of the power x in all powers o* 
where 1 = k = n. In a fixed power o*, the coefficient 


of x arises as the sum of all products 


h 
(~1) a a a ‘ 
a tp Ky 
where k. Sig. Ss Ay. 27 , hh, and 
Ky + kK. + + KD =n 


where the sum comprises all systems of indices Kye Kor 


, k, such that h ae 


IV 


(2) 


(No upper limit on h needs to be imposed.) 

At first sight, it seems difficult to program a 
complicated formula such as (1) for a pocket calcula- 
tor, due to the many systems of indices (kK, ,k5,---+7k,) 
that have to be kept track of. However, the following 


device, suggested to the author by E. Specker, makes 
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the programming possible. The systems of indices k = 


(kK) Ko --- rk) satisfying (2) are put into a one-to- 


one correspondence with the integers m satisfying 


ae S m < eae in the following way: Let 
2 n-Ll 
m= €), + 2c, +. 2 Ee, + a oe Cn-1 (3) 
be the binary representation of m oe = 0 or l, 3 = O, 
dey Bidder pe ZS Sa = 1). For each j such that = = 1, 
we put a mark at the point x = j of the real line. In 
addition, we put a mark at x = 0. These marks divide 


the interval [0, n-l] into subintervals. The indices 


kK. of the system (ki rkore ee rk 


the lengths of these subintervals. 


) corresponding to m are 


Accordingly, the coefficients bo may be calculated 


by the formula 


where the products I are defined by means of the bi- 


nary expansion of m, 


k k, are the gaps between the non- 


qi? iad eo © @ tf h 
Zero =e im: ¢3).2 


here k 


The binary expansion of m is calculated by the 


following algorithm. Let Py c= Mm, 
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Then 


O, otherwise. 


Given m, the product [I :=  ™) thus is calculated re- 
cursively as follows: Let yO) s= ly 
(3? es = 0, 
7 (att) 
Bie 956g POS 
here k := 3 - j' where 3' := max { i | i<j, E5 x 0 }. 
(n) 
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3. Flow Diagram 


4. 
R 


0 


Storage and Program 


Reciprocal Power Series 


CHS 
STO*6 
RCL 5 
x # 0 
GTO 07 
RCL 6 
STO+7 
RCL 2 
x # 0 
GTO Ol 
RCL 7 
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Locations 20 + 38 are reserved for the program for the 
computation of aye This program should assume k in R73 
at the end of program execution, ay must be in X re- 
gister. The registers Ro and R, are available for 
auxiliary storage. 

5. Operating Instructions 

Load the program, including the program to compute ay 


Unused locations must be filled up by NOP instructions. 
Move the operating switch to RUN. Select the mode of 
displaying numbers. If coefficient bo is desired, load 


data as follows: 


2 into Ry 
eas into R, 
0 into Ro 


Here the numbers 2” and ae must be genuine integers; 
the calculator will not deal correctly with a 2” that 


is computed by the 7 instruction. Press 


PRGM 
R/S 


to start computation. After some time (see examples) 
the calculator will display bo: Caution: bo is calcu- 


lated as a sum of 2” terms which may differ in sign 
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and size. Thus already for moderate n, bo may be af- 


fected by rounding error. 


6. Examples 


To calculate the coefficients bo in the expansion 


xX 2 
Tog Se ee ee + 
(The be are required for certain methods of nume- 
rical integration; see ENA, § 13.4.) The fore- 


going series is the reciprocal of 


x 2 3 


The coefficients 


are generated by means of the program 


20 RCL 3 24." 

21 1 : 

22 + NOP 
23 1/x 38 


Results: 


132 


Power Series 


bo (calculated) oo (exact value) 
-0.500000000 - = = -0.500000000 
-0.083333333 - 3 = -0.833333333 
-~0.041666667 - = = -0.041666667 
19 
-0.026388889 - =—— = -0.026388889 
720 
3 
-0.018750000 - Té0 = -0.018750000 
-0.014269180 a -0.014269180 
. 60480 : 


Computing 


time 


25:2 


Sec 


sec 


SEC 


sec 


sec 


SECC 


All calculated values are accurate in all digits 


given. 


To compute the Bernoulli numbers By 
x e. ae k 
oS) ae 
e - l k=0 ~ 
(see ACCA I, p. 13). - Setting 
em = 4 


Se SL. aa x a x? + 
xX 1 2 = 


the program will calculate the numbers 


defined by 
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The coefficients 


in ee ieee tae 
k (k+l)! 


may be calculated by the following program: 


20 1 28 GTO 23 
21 STO+3 29 RCL 3 
22 RCL 3 30 1/x 
23 1 31 
24 - 
25 x = 0 : NOP 
26 GTO 29 
27 STO*3 38 
Results: 
be (calculated) bo (exact values) Computing 
time 
-0.500000000 - - = -0.500000000 3 sec 
0.083333333 +3 = 0.083333333 8 sec 
-ll 
3 * 10 0 20 sec 
-~0.001388889 - as = -0.001388889 5A. see 
9 * 10 7? 0 129 sec 
0.000033069 L = 0.000033069 303 sec 
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POWER OF POWER SERIES 


1. Purpose 


Given a formal power series with leading coefficientl, 


2 
P=i1+ a,x + a,x + a,x ee ee (1) 


and given a real number a, to compute the coefficients 


a of the a-th power of P, 


2 
P = (1 + a,x + a,x ey setae) 


Il 
~~ 
+ 
s0)) 
m* 
+ 
03) 
mm 
ee 


2. Method 


By conventional wisdom, the purpose would be served by 
means of a recurrence relation derived from the fact 
that the series R := P” satisfies the formal differen- 


tial equation 
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(J. C. P. Miller formula; see ACCA I, Theorem 1.6c). 
This approach is not feasible for pocket calculators 
due to storage limitations. Instead, we compute the 
series P” directly from its definition as composition 


of the binomial series 


This yields 


hence 


where the sum comprises all systems of indices (kk 


paherpete 


pas 
h) satisfying the conditions (2) of the preceding 
algorithm. These systems are coded in terms of the bi- 
nary representation of the integers m satisfying gr 


< n 
=m < 2 , and we obtain 
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where ue has the same meaning as in the preceding al- 
gorithm, and where h is the number of the non-zero di- 


gits in the binary representation of m. 
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3. Flow Diagram 


Power Series 


. Storage and Program 


FG ay ae 

on ono 
m R 
00 

> O01 1 
02 STO 6 
03 STO-1 
04 STO-2 
05 RCL 1 
06 STO 5 
07 CLX 
08 STO 4 

= O09 CLX 
10 STO 3 

= del. 1 
12 STO+3 
13 RCL 5 
14 2 
15 ? 
16 INT 
17 STO 5 
18 LAST x 
19 FRAC 
20 x = 0 
21 GTO ll 
22 
23 
24 
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Locations 22 + 31 are available for the program to 


compute a This program should assume k in R at the 


‘ ; 
end of oe execution, ay must be in X a 
(If it is more convenient to compute lfay, this may be 
done; instruction 32 must then be changed to STO+6.) 
No registers are available for auxiliary storage, but 


stack may be used freely. 


5. Operating Instructions 


Load the program, including the program to compute ay: 
Unused program locations are to be filled with NOP in- 
structions. Move the operating switch to RUN. Select 


the mode of displaying numbers. Load exponent 


oO into Ro 


(a) 


and, if coefficient an is desired, 


2 into Ry 
2 into R 
0 into R 


(The powers of 2 should be loaded as integers, without 


using the y* instruction.) Press 


PRGM 
R/S 
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to start computation. The calculator will stop by dis- 


coe Caution: Computing time increases expo- 


n 
nentially with n. For large n, ae may be contamina- 


playing a 


ted by rounding error. 


6. Examples and Timing 


Power series reversion. Let y be given by a 


power series in x of the form 


Then the expansion of x in terms of y is 
Sy hey 
x = Y oY 3Y cee, 


where, by the Lagrange-Biirmann Theorem (see ACCA 


ry §, 1.9) 


Here res Rg denotes the coefficient of .- in 
-n 
the expansion of R=. Because R = xP, where P has 


the form (l), this coefficient in the foregoing 
(~n) 


notation equals a5 There follows 


setae 


n a ne. ar 3, 
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As an example, consider 


R=e -1l, 
where 
a = a 
k (k+1)! 
Because 
co n-l 
REU = Log(l1 + x) = ) het : 
n 
n=l 
we have in this case 
n-l 
a oe (-1) . Pa n) _ (2 1 
n n 


The following program is used to compute the ay: 


22 1 28 GTO 31 
23 STO+3 29 STO*3 
24 RCL 3 30 GTO 25 
25 1 31 RCL 3 
26 = 32 STO+6 
27 x = 0 


Being ignorant of the exact answer, the main pro- 


gram then computes the following values: 
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n ee Computing time 
2 -1.000000000 4 sec 
3 1.000000000 10 sec 
4 -1.000000000 23 sec 
5 0.999999999 58 sec 
6 -1.000000000 137 sec 
7 0.999999998 5.5 min 
8 -1.000000006 12 min 
9 1.000000008 
10 -0.999999858 
11 0.999999631 


Stirling's formula. The coefficients c, in the 


asymptotic expansion of the [I function, 


Cc Cc 


T (x) ~ /2T gene ae {l + = + -4 + .eet, KX 7%, 
x 
can be defined as 
1 
= 59: 7h \=qe-) 
“g 2 (Dg oq 2 ’ 
where 
P = , a ra = 1+ Ze + 2. + 20° 
z k 7 3 4 5 


(see ACCA II, § 11.6). The program calculates the 


following values: 


0.083333334 


0.001157408 


~0.000178759 


-0.000002234 
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0.083333334 


0.003472223 


-0.002681385 


-0.000234570 - 


2488320 


exact value 


1 _ 
Io = 0.083333333 
ie = 0.003472222 
288 : 
139 
- BS1840 7 -0 .002681327 
571 


-0.000229472 
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EXPONENTIATION OF POWER SERIES 


flees Purpose 


Given a formal power series with constant coefficient 


zero, 


to calculate the coefficients Ce of the series 


Q_ 2 
e~ = 1 + C yx + Cox + cee 


2. Method 


To use the differential equation satisfied by the se- 
ries S := ee is not feasible due to the lack of stora- 


ge facilities. However, from the fact that 


e~=1+ ) 
h= 


we have, using the notation of the two preceding pro- 


grams, 
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where the sum comprises all systems of indices (kk 


7k 


D4 


n) satisfying the conditions (2) of the program 


"Reciprocal Power Series". As before, these systems 


are coded in terms of the binary representation of the 


integers m satisfying ges =m < ae and we thus ob- 
tain 
gel 
Ee = 3° 2 
n hi om ’ 
n-l 
m=2 


where I has the same meaning as in the two preceding 


programs. 
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3. Flow Diagram 
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4. Storage and Program 


Ro Ry Ro R, Ry R- Re Ro 

a ai h p H |o| 
m R 2 
00 25 

> 01 1 26 
02 STO 6 27 
03 STO 4 28 
04 STO-1 29 
05 STO-2 30 
06 RCL 1 3) 
07 STO 5 a2 

> 08 CLX 33 
09 STO 3 34 

+ 10 1 35 
ll STO+3 36 STO*6 
12 RCL 5 37 RCL 4 
i Be 2 38 STO+6 
14 3 39 1 
15 INT 40 STO+4 
16 STO 5 41 RCL 5 
17 LAST x 42 x F 0 
18 FRAC 43 GTO 08 
19 x = 0 44 RCL 6 
20 GTO 10 45 STO+7 
21 46 RCL 2 
22 47 x # 0 
23 48 GTO 01 
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Locations 21 + 35 are to be used for the program to 


compute a This program should assume k in R at the 


k" aH 
end of program execution, ay must be in X register. 


(If it 1S more convenient to compute l/a,, this may be 
done; instruction 36 must then be changed to STO+6.) 


Register R, is available for auxiliary storage. Unused 


0 
locations are to be filled with NOP instructions. 


Bis Operating Instructions 


Load the program, including the program to compute aye 
Move the operating switch to RUN. Select the mode of 


displaying numbers. If coefficient Co is desired, load 


© 
~- 
) 
ct 
O 
ve) 


(The powers of 2 should be loaded as integers, without 


using the y* instruction.) Press 


PRGM 
R/S 


to start computation. The calculator will stop by dis- 
playing co: (Computing time increases exponentially 
with n.) Caution: For large n, if some ay are negative, 
cancellation may cause Cs to be contaminated by roun- 


ding error. 


6. 


Exponentiation of Power Series 


Examples and Timing 


Q = Log(l - x) 


We should obtain 


that is, Cy 


to compute a 


Results: 


mM F& W DN EE 


k 


See Dg 
k=1 . 
eX =ljl-x, 
atl Oar Cy = 0 for k > l. 


is as follows: 


21 RCL 3 
22 CHS 
23 GTO 36 
24 | 


: NOP 


36 STO=6 


Cc 
n 


-~1.000000000 
0.000000000 
0.000000000 

oe as 


0.000000000 


Computing time 


3 
6 
15 
38 
90 


sec 


sec 


sec 


sec 


sec 
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The program 
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6 0.000000000 3.0 min 


7 1 * 10 29 Sot aia 


The Bell numbers bo are defined by the expansion 


[o@) 
n nN 
e = y nt 


For 


we may calculate a, by the program 


k 


21 RCL 3 28 RCL 3 
22 1 29 STO=6 
23 = 30 GTO 37 
24 x = 0 31 
25 GTO 28 ; 
26 sTo*3 pear 
27 GTO 22 36 
Results: 
n Cc b =nic Computing time 
n n n 
1 1.000000000 1 3 sec 
2 1.000000000 2 8 sec 


3 0.833333333 18 sec 


1S) 


oOo On HD NN & 


oOo oOo Oo Oo OO © 
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-625000000 
433333333 
-281944445 
- 174007996 
~102678571 
-058275463 


15 

52 
203 
877 
4140 
21147 


46 
110 
260 

10 

23 

oh 


sec 
sec 
sec 
min 
min 


min 
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Part 5 
INTEGRA TION 
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NUMERICAL INTEGRATION WITH STEP REFINEMENT 


1. Purpose 


To evaluate the integral 


numerically for an arbitrary function f that can be 


evaluated by a program requiring no more than 16 in- 
b 

structions. The general integral f is reduced to the 
a 


above by the substitution x' = x a. 


2. Method 


We approximate I by a sequence of approximate values 
Ty, k = 0, 1, 2, ... . The approximation Ty is ob- 
tained by dividing [0, b] into 5% congruent subinter- 
vals and evaluating the integral on each subinterval 
by the midpoint formula. (The trapezoidal rule is 
avoided because integrands frequently have singulari- 


ties at the endpoints of the interval. Although harm- 
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less for the existence of the integral, these may 


cause the program to fail; see example .) Writing 


A ee DO Sy : 
we have 
eal ; 
I = hy z £((m + 5) h,) 
m=0 
For any continuous f, 
lim I. = I 
ko 


If f is sufficiently smooth, the convergence of the 
sequence {IJ may be sped up by the Romberg algorithm. 
To this end the program saves the five currently most 


recent values of I, in locations in which they can be 


k 
used as input for the Romberg acceleration program 


(see the following program) without external storage. 
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3. Flow Diagram 


compute £(x) 
s := s + f(x) 


x +h 


Display 


push down I 


1 
h := 5h 


k 
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4. Storage and Program 


m3 PA 
I, Ty 
(R/S) 


DOD MN WwW FP NYO WH FF WD 
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The program to compute f(x) should be in locations 33 


through 49, assuming x in R Only the stack may be 


6° 
used for temporary storage. The last instruction must 
be GTO 09. At this point, f(x) must be in the X regi- 


ster. 


5x Operating Instructions 


Load the main program. Load the program for computing 
f into locations 33 through 49; the last instruction 
must be GTO 09. Switch to RUN. Select a mode of dis- 


playing numbers, for instance 
FIX 8 


If f£ involves trigonometric functions with argument in 


radians, press 


To start computation, press 


PRGM 
R/S 


The calculator will pause briefly while displaying 


each I,. If k 20, and if 


R/S 


160 Integration 


is pressed during pause, the values ee will be found 


in Rayan! m= "0..-- deg: 2k. iy Ae 


6. Examples and Timing 


bs Sei AEs): = A, . Program to compute f: 
1+ x 
33 4 
34 RCL 6 
35 * 
36 1 
37 + 
38 + 


39 GTO 09 


The following values of I, are obtained: 


k 
(k) 
k IV I, 
0 3.20000000 | 3.14191817 
1 3.16235294 3.14159265 
2 3.14680052 Romberg, 3.14159264 
3 3.14289473 3.14159266 
4 3.14191817 | 3.14159264 


Time required 35 sec. Subjecting the values T. 


to the Romberg algorithm produces the exact value 


T = 3.14159265 with an error of one unit in the 
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last digit due to rounding effects. 


=) HOD Net) . Program to compute f: 


X 
33 1 
34 RCL 6 
35 + 
36 ln 
37 RCL 6 
38 + 
39 GTO 09 
Values obtained: 
(k) 
k Ty I, 
0 0.81093022 ' 0.82241711 
1 0.81936429 0.82246693 
2 0.82167416 Romberg. 0.82246702 
3 0.82226766 0.82246702 
4 0.82241711 0.82246702 
ne 
Exact value: ID 0.82246703. Time required to 
generate Ig: gttee i I, about 40 sec. 
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ROMBERG ALGORITHM 


1. Purpose 


To accelerate the convergence of a sequence iT} to 


its limit Ags where 


a T(2 "h) , m=0, 1, ...,4, 


and the function T(h) has an asymptotic expansion of 


the form 


2 4 6 
y 
Ph) Wee tae eh ean, ecg J Oe 


This problem poses itself in connection with numerical 
integration [T (h) = approximate value of integral cal- 
culated with step h]; in addition, it occurs in nume- 
rical differentiation, and in the computation of cer- 


tain functions and constants. See examples below. 


2. Method 


The Romberg algorithm (see ENA, § 12.4; ACCA II, 
§ 11.12). With the data 
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TO t(2 “h) , m=0, 1, ..., 4 


we generate a triangular array of numbers T an by the 


formulas 


T ee Ie my ly eee Se A 


If m were allowed to tend to infinity, one would have 
(loc. cit.) for each fixed n 


T sa +0(4™) 


Thus Tag may be regarded a much more accurate approxi- 


mation to Ao than TAQ (mathematical error bounds are 


available). To check convergence, the program is ar- 


ranged so that all values T (n = O, 1, ... , 4) are 


An 
available at end of computation. 
For reasons of numerical stability as well as eco- 


nomy, formula (1) is evaluated as 
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3. Flow Diagram 


Because an address modification is not available on 
the HP-25, the operation (2) is always carried out at 
the same place. This is made possible by a cyclic per- 


mutation (here called "rotation") of the Ta 


Yes 


Rotate 


:= m+tl 


Display T 


Stop 


4. 


+ 03 


+ 05 


+ 23 


ys) 
Q 
tc 
wKe IN WD NI 


RCL l 
+ 
STO 0 
PAUSE 
RCL 0 
STO 5 


Storage and Program 


(R/S) 


Romberg Algorithm 


- O W FPF NH W FF NY CO FF 


03 


00 
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5. Operating Instructions 


Load the program, and switch to RUN. Choose the de- 


sired mode of displaying numbers, for instance 
SCI 8 


to get floating eight-digit display. Load data ae = 
T(2 “h) into Rw! m= 0, l, ..., 4. (If the data are 
produced by the program "Numerical Integration with 


Step Refinement", they are already there.) When 


PRGM 
R/S 


is pressed, the calculator will generate the array Des 
column by column, pausing briefly after computing each 
entry. (For an indeterminate display of each Tait 

change instruction 22 to R/S and press R/S after each 


display.) The calculator will stop by displaying T 


44° 
Ta4 is also in Roi generally, the elements in the last 
row of the array, Ty n’ are found in the registers 
a 
Rant he Oop, ee - ete. op, “Ae 


6. Examples and Timing 


Numerical evaluation of 
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by the preceding program produces the approximate 


values 


0.001953125 
0.037544250 
0.078839094 
0.094288322 
0.098544474 


The present program accelerates the convergence 
of these values to produce the following values 
Pa of in the last row of the Romberg array (listed 
vertically here for reasons of space): 


0.098544474 
0.099963191 
0.099998199 
0.099999858 
0.100000000 


It is seen that the last value is accurate to all 


places given. Computing time about 44 sec. 


[2] Here we compute Euler's constant, 


(3 aie Se = @ 


1 
5 ? Gedue, an Log n) , 


Wir 


n> 


168 
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directly from its definition, using Romberg acce- 


leration. If for n= 1, 2, ... we let 

1 1 1 11 
s(n) :=1l+5 4+ 35+ .-- +57 +577 bog nn, 
then we evidently also have y = lim s(n). The 


factor - in front of = is motivated by the fact 


that for suitable a 


(see ACCA II, § 11.11). Thus the numbers 
k 
G. t= s(2) , Ko Oy. Le 27: «ae | 
satisfy the conditions for applying the Romberg 


algorithm (see ACCA II, § 11.12). They are con- 


veniently generated by the recurrence relation 


Sea. 


k = 0, 1, 2, ... . It is an easy exercise to 
write a program that computes o. for i= 0, l, 


, 4 and stores it in Ry. We obtain 
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-500000000 
-556852820 
.5972038973 
-575915602 
-576890272 


(o> es © ee © ee ee) 


The program "Romberg Algorithm" accelerates these 


values as follows: 


0.576890272 
0.577215162 
0.577215652 
0.577215663 
0.577215665 


Again the last value is accurate to all places. 


Computing time for the acceleration about 44 sec. 


For further examples on Romberg's algorithm see 
the programs "Numerical Integration with Step 
Refinement", "Plana Summation Formula", "Log- 


Arcsine Algorithm". 
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PLANA SUMMATION FORMULA 


1. Purpose 


To evaluate the integral 


y 
ec TY a 


L Se Pee a 

0 
by means of numerical integration. This integral oc- 
curs, among other places, in the Plana summation for- 


mula (ACCA I, § 4.9), 


5 £(n) = =£(0) “PEG ax 4.7 ah ay 
- 0 0 e - l 


where 


g(y) := - 2 Im fliy) 


This formula holds, roughly, for functions f that are 
analytic for Re z = 0 and grow less than e2T ly! for 
y := Imzwssto~ , It may often be used to sum slowly 
converging series. The program works for any function 
g whose evaluation requires no more than 18 instruc- 


tions. 
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2. Method 


The integral I is approximated by the sums 


(00 ) gly,) 
I. = ah , m= 0, Ly, 2, ; 
k=0 Yk 
e - l 
-m h ; 
where h=h, - 2, y, = = + kh. The summation is 
0 Ty k 2 
stopped when e - 1 > 1l/e. The constants Ao and l/e 


are input. The program displays each L. after it has 
been computed. The convergence of the sequence {1} 
may be accelerated by the Romberg algorithm (see the 


preceding algorithm). 
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3. Flow Diagram 


4. Storage and Program 


ANV 


* NO 


RCL 6 


Plana Summation Formula 
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Note: The function g is evaluated by instructions 08 
through 25, presuming its argument y in Ree gly) 
should be in X register before executing instruction 
26. In above program, g is taken from example be- 


low. 


5. Operating Instructions 


After loading the program, switch to RUN. Indicate the 


mode of displaying numbers desired, for instance, 
SCI 8 
If function g involves trigonometric functions with 


arguments in radians (such as in examples (2 and 


below), press 


RAD 
Load Ay (such as 0.5) and press 
STO 0 
Load 1/e (such as 10" = 10) and press 
STO 1 


When 
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PGRM 
R/S 


is pressed, the calculator will compute and display 
I.. Pressing R/S again yields I ai 


0 LW arotee 
The numbers I 


0! Ij are not stored. To subject 
them to the Romberg algorithm, they should be copied 


and used as input for the program Romberg Algorithm. 


6. Examples and Timing 


[1] Euler's constant y (see preceding algorithm) may 


also be represented as 


veges E> a wy 
0. Bay ae ea A 


(ACCA I, p. 275). The function 


2 
hy) 42S 
l+t+oy 
is generated by 
08 RCL 6 
09 2 
10 = 
11 RCL 6 


12 x 
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13 

14 + 

15 + 

16 GTO 26 


The following values I, result by choosing ve := 
0.54 1/6 2 aaa 


.6297262 * 10 
.4582004 * 10 
.6562828 * 10 
.7052792 * 10. 
7.7174967 * 10. 


~“Jlo o~! wn OD 


Total computing time: approx. 10 min. Subjecting 
the foregoing values to the Romberg algorithm 
yields the following last line of the Romberg 


scheme: 
7.7174967 7.7215692 7.7215663 7.7215664 =7.7215664 


(all * o>). We conclude that I = 7.7215664 * 
10>, and thus 


y = 0.577215664 , 


which is correct to within a rounding error of 


i he aes 
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Application of the Plana formula to f(z) := 
(l + z) ° (principal value) shows that the Rie- 
mann zeta function (ACCA II, § 10.8), 


ECs): 2= b aay Le 


may be represented as follows: 


[0 


1 1 
ELS): = She — - + f g(s;3 y) dy , 
2 S 1 0 er TY _] 
where 
_ 2 sin(s Arctg y) 
g(s; y) im 2 s/2 
(oak ye) 


This representation holds for all s # 1. We load 


Ss into R, and use the following program to com- 


3 
pute g: 
08 RCL 6 
09 can 
10 RCL 3 
11 * 
12 sin 
13 2 
14 * 
15 RCL 6 


16 x 
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17 1 
18 + 
19 RCL 3 
20 y 
21 vx 
22 = 


For s := 1.1, Ao := 0.5, lfe := 107° we obtain 
712387260: * 10 °°) ( 8.4403697 * 10° 
8.1548623 * 10° 8.4448497 * 10-7 
8.3730175 * 10° Romberg. < 8.4448464 # 16° * 
8.4269295 * 107° 8.4448464 * 107° 
8.4403697 * 107° 8.4448464 * 10° 


There follows 7(1.1) = 10.584448464. The program 
may be tested by evaluating ¢(s) at the "trivial" 


zeros s = -2, -4, 
To evaluate the very slowly converging series 


ee ee 
2 
3 n{Log n} 


if 


S := 2 
n= 


we use the Plana formula where 


1 


£(z) := ee ee 
(3 + z){Log(3 + z)} 
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This requires the evaluation of 
g(y) := - 2 Im f(iy) , 


which may be written as 


_ 2 sin(a + 28) 
rR 


where 


a Ly eS e% , Log r + ia = R ae 

The program is given in the foregoing main pro- 
gram; it offers interesting applications of the 
> P operation. The final multiplication by 2 is 


omitted for lack of space. The following values 


result with hy S=-102. 5. Lee = eae 
~2 y ( = 

0.8983867 * 10 1.0359638 * 10 
1.0025014 * 10°? 1.0364923 * 107 
1.0280286 * 10° Romberg. 1.0364921 * 10° 
1.0343783 * 10° 1.0364921 * 10 7 
1.0359638 * 10° 1.0364921 * 10° 
We conclude 

= 1 ee es ae ee * * —2 
S = 5 et ears ey 1.0364921 * 10 


3(Log 3) 


= 1.0690583 . 
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The total time required for this computation is 


about 20 min. 
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DIFFERENTIAL EQUATION OF FIRST ORDER: 
TRAPEZOIDAL METHOD 


ie Purpose 


To integrate numerically a general first order equa- 


tion 
yr = E(x ey) (1) 
Subject to the initial condition 


y (X_) TG 


2. Method 


The trapezoidal rule. Choosing an integration step h, 


we determine approximate values of the values y(x_) 
Yn n 


of the exact solution y(x) at the points x t= Xp + nh 
(n= 1, 2, ...) by the formula 
~y = Sle(x vy) + £( ) (2) 
Yn+1 Yn 2 n’*n *nt+1'¥n4+1 : 
The implicit equation for y that arises at each 


n+l 
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step is solved by iteration, using Euler's method as a 
predictor (see ENA, § 14.7). Iteration is stopped when 
two successive iterates differ by less than ¢«. Eleven 
storage locations are available to evaluate the ex- 

De (ei). 


2 
The order of the method defined by (2) is 2, inde- 


pression 


pendently of the number of iterations performed. That 
is, as h > 0, the approximate values ve at a fixed va- 
lue of x tend to the exact value y (x) with an error 
of O(h*). If « is small (say, «= 107°) and the equa- 
tion (2) thus is satisfied to full machine accuracy, 
the convergence of ve to y(x,) may be sped up by the 
Romberg algorithm (see the example). However, a large 
number of iterations may then be necessary at each 
step. The number of iterations may be reduced by gi- 
ving € larger values. For very large values of ¢€ (Say, 
Ee = to the corrector will be applied only once, 
and the method becomes identical with the simplified 


Runge-Kutta method, 
= _h f(x ) + £(x +hf£ ( )) (3) 
Yn+1 Yn 2 n'*n n+1'%n Yn ; 
The order of the method is still 2, but Romberg's ho 


extrapolation is no longer feasible. 


3. Flow Diagram 


We write x := 2 ae y* denotes the current ar- 
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gument y in £(x,y), yo is the corrected value. 


Display x 


Display y 


A "flag" k is used to distinguish between the predic-— 


tor (k = 0) and the corrector cycles (k = 1). 
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4. Storage and Program 


+ 05 


Program to 


compute 


£(x,y*) 


RCL 0 
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5. Operating Instructions 


Load the main program. Load the program for computing 
f(x,y*) into program locations 05 + 15, assuming x in 


R, and y* in R This program must place f(x,y*) into 


- (As an ied panda possibility, one may use locations 
05 + 20 for a program that places A := 2 (x,y*) into 
Re -) R, is available to store a constant required for 
the evaluation of f. 

Move the operating switch to RUN. Select the de- 
sired mode of presenting results, for instance by 


pressing SCI 6 to get floating six-digit display. Load 


integration step h into Ro 
starting value Xo into Ry 
starting value Yo into R, 


tolerance ¢« for stopping 


, : into R 
inner iteration 3 


When 


PRGM 
R/S 


is pressed, the calculator will start computing. After 
completing inner iteration, it will briefly display xy 
and stop by displaying Yi° Pressing 


R/S 
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will cause the display of x etc. The process may 


9' Yo" 
be repeated until the desired range of x has been co- 
vered. 

To check the accuracy of the results, the computa- 
_ oe 
til the results have converged to the desired accura- 


tion should be repeated with the steps 


cy. If € is small, the convergence may be sped up by 


the Romberg algorithm. 


6. Example and Timing 


The exact equation of the mathematical pendulum of 


length 2&, 


(g := gravitational constant = 9.81 ms“) may under 


the initial condition y'(0) = 0 be integrated to give 


y' = - /=3 Ycos y - cos Ne 


where Vg y(0). We solve this equation under the ini- 


tial condition 
= y(0) = = = 1.047198 
Yo ~ ¥ ~ 3 


Here 
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2a oS Ss OD mana @ as OS 
——— 5E (x,y) = h x0 cos y 0.5 


The following program is used to compute A: 


05 RCL 7 
06 cos 
07 

08 5 
09 - 
10 Vx 
11 RCL S5 
12 * 
13 RCL O 
14 * 
15 CHS 
16 STO 6 
17 NOP 
18 NOP 
19 NOP 
20 NOP 


The constant / <2 is assumed in Re Since the argument 
of cos is assumed in radians, RAD must be pressed be- 
fore starting computation. 

The following values result if 8 = 1, beginning 


with h = 0.1 and successively halving the step: 


Integration 


h=0.1 h=0.05 h=0.025 h=0.0125 
0.0 1.047198 1.047198 1.047198 1.047198 
O.1 1.005241 1.004979 1.004913 1.004896 
0.2 0.881649 0.880558 0.880284 0.880215 
0.3 0.683836 0.681290 0.680648 0.680488 
0.4 0.425723 0.421160 0.420008 0.419719 
0.5 0.128475 0.121689 0.119976 0.119546 
0.6 -0.180863 -0.189433 =0i2.4, 95159 5 -0.192137 


Applying Romberg to the values at x = 0.5 yields 


0.128475 

0.121689 0.119427 

0.119976 0.119405 0.119404 

0.119546 0.119402 0.119402 0.119402 
For h = 0.1 and an iteration tolerance of ¢€ = 


co. the time per integration step varies between 15 


and 60 sec, due to the considerable number of itera- 
tions required. By setting € = Lo the program is 
modified to correct only once, possibly at the expense 
of stability. The time per integration step is then 


reduced to 6 sec. 
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AUTONOMOUS DIFFERENTIAL EQUATION OF SECOND ORDER, 
FIRST DERIVATIVE ABSENT 


bg Purpose 


To integrate numerically a differential equation of 


the special form 


ye EY). (1) 
subject to the initial conditions 
y (Xp) YQ ’ Y (xX) = 29 
Equations of the special form (1) occur, among other 


places, in problems of non-linear vibrations without 


damping. 


2. Method 


A simplified version of the Runge-Kutta method (see 


ENA, § 14.5). Let h be the integration step, and let 


0 
of y' (x) are computed by the formulas 


= + , 
x xX nh. Approximate values ve of y(x) and Z 


190 


3 


Integration 


Oy. by 25 


x is fixed is 


Flow Diagram 


We write x := x _, 


current argument in f; 


nN 


(2) 


h 
= Z + hily, + 525) ; 


The error in y.andz ash-> 0 and 
5 n n 
O:CBi +) 


Viger ee Zi y* denotes the 
A 


= ee Cys) 
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Display x 
Display y 
Stop 


192 


4. Storage and Program 
Ry Ry Ro R 
h x Y Z 

00 

Ol CLX 

02 STO 4 

03 RCL 2 

04 STO 7 
+ 05 ) 

06 

07 

08 

09 

10 Program to 

11 

12 compute 

13 

14 f(y*) 

15 

16 

17 

18 

19 

20 RCL 0 

Zi - 

22 STO 6 

23 RCL 4 


Integration 


24 x ~# 0 
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5. Operating Instructions 


Load the main program. Load the program for calcula- 
ting f(y*) into program locations 05 + 19, assuming 
y* in R.. Program must put f(y*) into X. Re is availa- 
ble to store a constant needed for the calculation of 
Fe 

Move the operating switch to RUN. Select a mode of 
displaying numbers. If £ contains trigonometric func- 


tions with arguments expressed in radians, press 


RAD 
Load data as follows: 
h into Ro 
Xo into Ry 
Yo into R, 
Zo into R, 


Put auxiliary constants into R. Pressing 


will start calculation. The calculator will briefly 
display xy and stop by displaying Y.° To exhibit Za 
press 


RCL 3 
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To continue computation, 


This will cause display of x 


6. 


Integration 


Example and Timing 


R/S 


2 


press 


and Yo: etc. 


We integrate the exact equation of the mathematical 


pendulum, 


where g/ 


ee os ae 


y (0) 


= ' Siny , 


ut 
3 a 


y' (0) 


= 0 


under the initial conditions 


(compare the program "Trapezoidal Rule"). The results 


for various values of h are as follows: 


Oo 2 a2 Oo Oo © "© 
DN WM F WD NY FH CO 


oOo Oo OF OO FF FF 


.047198 
.004719 
.878363 
-675164 
-408987 
. 102068 
.216063 


Oo 2 © 2&2 — 


h=0.05 


.047198 
.004785 
-879600 
-678906 
-416708 
-114850 
-198319 


Oo “Oo Oe SO fe oF 


h=0.025 


.047198 
.004859 
-880033 
-680034 
-418874 
-118251 
-193810 


h= 


CO O. “OO. 


0.0125 


.047198 
-004885 
-880156 
- 680340 
-419443 
~119124 
- 192679 
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Applying h?-extrapolation to the values at x = 0.5 


yields 


.102068 
~114850 
~118251 
~119124 


oOo CO GO O&O 


0.119110 
0.119384 
0.119415 


The theoretical hypotheses for continuing Romberg ex- 


trapolation beyond the second column of the Romberg 


scheme are not satisfied here. 


Computing time per integration step (including 


display of x) approximately 5 sec, independently of 


choice of h. 
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LINEAR SECOND ORDER DIFFERENTIAL EQUATION 


1. Purpose 


To integrate numerically the general linear second 


order equation, 


y+ p(x)y” + oq (a)jy =O 5 (1) 
Subject to initial conditions 
Y(X%)) = Yo » y'(X%)) = Yq - (2) 


2. Method 


Finite differences. Denoting by Y, 2n approximate va- 
lue of the exact solution y(x) at the point x, t= Xp + 
nh, where h is the integration step, we approximate 


Y = oy 
; n+l] n-1l 
Y () ny 2h 


and 
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Introducing these approximations in the differential 


equation and writing 


B. 82 Pie) 4 @. a= a): y 
there results 
ey - 2y +y yee aDiy a )+qy = 0 
2 n+l n n-1l 2h 7 n+l n-1 n-n ’ 
which can be solved for Vogt to yield 
_ 1 2 h 
Yn+1 h 2 in Fae ae ot 2 A ae! 
Lr Sp 
2° n 
(3) 
_ 1 2 
— 2+ pp. ‘4 an ed BP Masa? 


This is used to calculate the approximations Y, recur 
Sively. The starting value Yo = y(Xq) is given. To ob- 


tain Y,, we use Taylor's expansion in the form 
2 


= ' laiaeeee 
¥y = Yo + BY9 * 2-¥o 


Here Yo and Yo = y' (X%q) are given, and Yo can be cal- 


culated from the differential equation. There results 


2 
Z © BL 
YMG RVG oe Povo Sava. oy 


This calculation must be performed manually. 
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3. Flow Diagram 


Compute y by (3) 


n+l 


Display x 


Display Vo 


Stop 
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4. Storage and Program 


Ro 4 PD ee Ry Rs Re aa 

h Xy Yo Yi temp temp 
*n Yn-1 Yn 
00 25 
O1 | 26 RCL O 
02 27 x? 
03 Program to 28 us 
ae compute oe # 
05 30 3 
06 P(x,) 31 CHS 
07 32 4 
08 33 + 
09 RCL 0 34 RCL 3 
10 ss 35 STO 2 
ll STO 4 36 i 
12 2 37 RCL 5 
13 RCL 4 38 = 
14 - 39 RCL 4 
15 RCL 2 40 2 
16 * 41 + 
17 STO 5 42 = 
18 | 43 RCL 
19 44 RCL 1 
20 Program to 45 + 
21 eonpeee 46 STO l 
22 47 PAUSE 
23 q(x) 48 R + 
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Sis Operating Instructions 


Load the program, including the programs to compute 
p(x) (locations 01 through 08) and q(x) (locations 18 


through 25), assuming x in R The registers R, and R 


are available to store sone wee needed in ee 
lation of p and q. 

Move the operating switch to RUN. Select the mode 
of displaying numbers, and choose correct mode for ar- 
guments of trigonometric functions (ordinarily by 
pressing RAD) if such functions occur in p or q. Load 


data as follows: 


h into Ro 
xy ( =X) +h) into Ry 
Yo into R, 


Calculate Yy from (4), and load it into R,- Pressing 
PRGM 
R/S 


starts computation. The calculator briefly displays Xo 


and stops while displaying Y>° Computation is conti- 


nued by pressing 
R/S 


which results in display of x Yor etc. Continue in 


a? 
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this fashion until the desired range of x values has 
been covered. 


As always with differential equation problems, it 


is a good practice to repeat calculation with h := sh 
in order to check convergence of values ae at fixed 
values of x. Romberg acceleration may be used at least 


once. 


6. Example and Timing 


We obtain values of the Bessel function J, (x) by inte- 


grating the differential equation satisfied by it, 


1 
A a oe ars ee a 


using 


NIF 


Yo = J, (0) = 0 , v5. = J, (0) = 


No harm is done by the singularities of the coeffici- 
ent functions p and gq at x = 0, because the program 
uses these functions only from x = xy onward. Equation 
(4), however, cannot be used to find Yi°- From the po- 
wer series defining J, (x) we get to the same order of 
accuracy 


oe? 
Y1 ~ 2 
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The programs to calculate p and gq are 


O01 RCL 1 
02 1/x 
and 
18 1 
19 RCL 1 
20 ee 
21 1/x 
22 - 


The remaining locations are to be filled up with NOP 


instructions. (As an alternate possibility, the whole 
program could be compressed, which is not problematic 
here because there are no internal references. In this 
case, the instructions R/S and GTO 00 should be in- 
serted at the end.) 

Values obtained using several different steps h 


are as follows: 


y J, (x) 

(exact 

x h=0.1 h=0.05 h=0.025 values) 
0.0 0.000000 0.000000 0.000000 0.000000 
0.1 0.050000 0.049958 0.049944 0.049938 
0.2 0.099667 0.099553 0.099517 0.099501 
0.3 0.148603 0.148406 0.148345 0.148319 
0.4 0.196436 0.196150 0.196063 0.196027 


=e W W WW W 


Computing time per integration step approx. 


° e e ° e 
oOo OO WO I BD 


-242806 


0.287368 


.329790 


.095396 


0.053631 


.012490 
.027698 
.066613 


Linear Second Order Equation 


242429 


-286899 
- 329230 


.095459 


0.053789 


.012740 
-027360 
.066193 


0 


~-242315 
-286758 
- 329063 


-095466 


0.053824 


-012801 
027274 
.066082 


203 


0.242268 


-286701 
- 328996 


-095466 


0.053834 


~012821 
-027244 
.066043 


4 sec. 


Part 6 


SPECIAL CONSTANTS 
ANAND FUNCTIONS 
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LOG-ARCSINE ALGORITHM 


Liss Purpose 


To compute the functions Log x and Arcsin x for arbi- 
trary x > 0 and for arbitrary x ¢€ [-1, Te respective- 
ly, by an algorithm requiring only the evaluation of 


square roots. 


2. Method 


Numerical differentiation (see ACCA II, § 11.12). 


Clearly, for x > 0, 


Log x = ae xh 
oh h=0 
Therefore, if 
Q(h) := s(x" - oy 
and 
n 
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then 
lim he = Log x 
n> 
Similarly, if x € [-1, 1| and a := Arcsin x, then 
Arcsin x = ae sin(ha) 
oh h=0 
Thus by putting 
— sin(ha) 
a(h) := ner eer od 
and 
a = a2} , n= 0, 1, 2, ; 
we have 
lim a = Arcsin x 
noo 
The sequences s_ = £ ands = a_ both satisfy the 
n n n n 
same recurrence relation 
2s 
— nN — 
“ol. “ae. + 6 a : 
n-l 


if 
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irae 1 1 1 
ai = q \X a =) r So ie 5 (x = z ’ (2) 
X 
then s, = 2 , and if 
Si = xvl - ie r Sp t= xX, (3) 
then Ss, = 4n: In both cases, the convergence of the 


sequence {s} can be sped up by the Romberg algorithm. 
The foregoing algorithm fails for the trivial values 
Log 1 = 0 and Arcsin 0 = O, but it is very stable 


otherwise. 


3. Flow Diagram 


The program computes {2} if the index k is set equal 
to 0, and it computes {a} if k # 0. In either case it 
saves the five most recent elements of the sequence 


{s_}. 
nN 
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Pause and display So 


Log-Arcsine Algorithm 211 


4. Storage and Program 


Ry Ry ie R, Ry Re Re Ry 
Sy S_3 S_» Sy So x 
00 25 1/x 
01 STO 6 26 - 
02 RCL 27 4 
03 x = 0 28 es 
04 GTO 15 29 STO 3 
05 1 + 30 RCL 1 
06 RCL 6 31 sto 0 
07 STO 4 32 RCL 2 
08 33 «sto 1 
09 - 34 RCL 3 
10 VX 35 STO 2 
11 RCL 6 36 RCL 4 
12 * 37 sto 3 

13 sTo 3 38 2 
14 GTO 30 39 * 

+ 15 RCL 6 40 RCL 3 
16 ENTER 41 RCL 2 
17 1/x 42 + 
18 - 43 + 
19 2 44 Vx 
20 + 45 RCL 3 
21 sto 4 46 * 
22 RCL 6 47 sT0 4 
23 2 48 PAUSE 
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5. Operating Instructions 


Load the program, and switch to RUN. Select the mode 


of displaying numbers, such as 
SCI 8 


to get the floating eight-digit display. Load k into 
Roy x into X (x #4 1if k = 0 and x # 0 if k # 0). 


Pressing 


PRGM 
R/S 


will start computation. The calculator pauses briefly 
while displaying each Su: No convergence test is pro- 
vided; convergence must be determined by inspection. 


Press 
R/S 


during pause to recover five most recent values. These 


will be in Ro through Rye 


the Romberg algorithm. The most recent value is also 


the correct locations for 


displayed in the X register. 


6. Examples and Timing 


k := 0, x = 10 yields with FIX 9 display 
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Log 10 = 2.302585094 in 16 iterations. 


Time required about 30 sec. There is an error of 
one unit in the last digit. By using SCI 8 dis- 
play we get 


Log 100 = 4.6051702 in 13 iterations, 
Log 1000 = 6.9077552 in 16 iterations. 


w 
7 
= 
x 
fT 


e yields the elements s 


ae) 0’ ’ 
S, of the sequence is} as follows: 

1.813430203 ' 1.002606203 
1.175201194 Romberg 0.999991848 
1.042190612 ier akg 1.000000050 

yields 

1.010449268 | —————_———-> | 1.000000001 
1.002606203 1.000000002 


Clearly, a FIX 8 display would have given the 


correct answer Log e = 1.00000000. 
k = 1, x := 1 yields with FIX 9 display 
Arcsin 1 = 1.570796319 in 15 iterations, 


which differs from the exact value 1/2 = 


P570796307 bY * Or. 
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THE GAMMA FUNCTION 


is Purpose 


To compute the values of the gamma function, 


r(x) 


ll 
k- 
s 


=fe t dt (x > 0) 
0 


for arbitrary real values of x with a relative error 
- es 

2. Method 

Stirling's formula (see ACCA II, § 8.5). We select n 


such that Stirling's formula evaluates [(x+n) with the 


required low relative error and then obtain 


(see ACCA II, § 8.4). Stirling's formula is usually 
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given in the form 


ris) 2 soe o (Xx - ¥2)Log x - x + J(x) (1) 


where J(x) is the Binet function, 


Log ee ee, 


-2 
t + x l-e i 


x > 0. It is known that J(x) possesses an enveloping 
asymptotic expansion aS x > ©, 


1 1 1 
1) a ee ee eke (2) 


360 ee 1260 ae 


The direct implementation of Stirling's formula, 
approximating J(x) by the terms of the asymptotic se- 
ries shown above, would produce the required accuracy 
for xtn > 5, but the resulting program would be too 
long for the HP-25. We therefore modify the formula, 


as follows. Writing (1) in the form 


ee (20 eX (Log x - £(x)) . (3) 


we see that it is necessary to evaluate the function 


fe). S40, = TU ; (4) 


which within the prescribed tolerance may be approxi- 


mated by the polynomial in hse. 


216 Special Functions 


1 1 1 
g(x) r= 1 - 5+ - (5) 
12 x 360 x 1260 x 


(see also Table l). 

To evaluate g(x) directly by means of (5), using 
any of the available methods for computing the values 
of a polynomial, would require the storage of three 
constants with a total of 9 digits, occupying 9 pro- 
gram memory locations. The resulting program would 
again exceed the capacity of the HP-25. 

It turns out that a shorter program is obtained by 
converting the function g(x) into a continued fraction 
(i.e., by representing it by its diagonal Padé appro- 
ximant), using the gqd algorithm (see ACCA, § 12.4). 
The following qd table results from the four coeffi- 


cients of the function g: 


(n) (n) (n) 
an 4a Sal To 
a 
ae 12 = 
12 lt 20 53 
1 30 58 315 
360 210 
_ 2 
_ i. 7 
1260 


This yields the continued fraction 
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as 
2 


© 


5 


WWJ 
1} Go 


which agrees with g(x) up to an error of at most 


O(x ~) as x + », In the above fraction we may replace 


the unwieldy = by = = - without affecting the ac- 


curacy of the final result. (Fe would arise exactly if 


7. : ae a ee 
the coefficient 1260 in g(x) were replaced by >1600 av 


55°) By algebraic manipulation the resulting appro- 


Ximation to £(x), 


may be simplified as follows: 


1 
Le 


h(x) = 7 


SKCAX: + ge eee se) 
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This, finally, is the approximation to f(x) used in 
our program. Only 4 one-digit constants need to be 
stored. The numerical values given in Table 1 show 
that the error committed in approximating g(x) by 
h(x) for x 2 5 again lies in the permitted range. If 
x < 5, Our program thus will determine the integer n 


< 
mentioned initially such that 5 = x+n < 6. 


X f (x) g(x) h(x) 

1 0.918938533 0.918650794 0.918566776 
2 0.979329652 0.979327877 0.979327701 
3 0.990774025 0.990773946 0.990773914 
A 0.994802332 0.994802324 0.994802324 
5 0.996671061 0.996671060 0.996671061 
6 0.997687312 0.997687312 0.997687312 
7 0.998300470 0.998300470 0.998300470 
8 0.998698592 0.998698592 0.998698592 

Table l 


3. Flow Diagram 


We let y := xtn, p := (x), - 
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No 


220 


Special Functions 


4. Storage and Program 


R 


0 


4 


00 
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5:4 Operating Instructions 


Load the program; operating switch to RUN. Select a 


mode of displaying numbers, ordinarily by pressing 


SCI 8 


to get 8 digit floating display. Load x into X regi- 


ster. Pressing 


PRGM 
R/S 


will cause the value of [(x) to be displayed. Values 
are generally correct to 8 significant digits; for 
large x errors may be slightly larger due to inaccu- 
racies in exponential function. 

To get [(x) for new value of x, load new x into X. 


Pressing 


R/S 


will yield new value. 
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6. Examples and Timing 


yields r(x) = 1.0000000 
1.0000000 
2.0000000 
6.0000000 
2.4000000 * 10 
1.,2000000 * 10 
7.1999999 * 10 
5.0400000 * 10 
4.0319999 * 10 
3 


-6287999 * 10 


Oo OO WO ND OHO FPF W DY FF 
ON £& W NY NN F 


oe 


These values provide a genuine check, because the 
relation [(n+l) = n! is not used in the program. 


Computing time for each value about 5 sec. 


[2] x = 0.5 yields r (5) = 1.7724538, which agrees 


with the exact value, /7. 


[3 | x = 0 correctly yields "Error", because [(x) has 
a pole at x = 0. The same holds for all negative 


integers. 


[4] x = -2.37 yields f(-2.37) = 1.2183595. 
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INCOMPLETE GAMMA FUNCTION 


l. Purpose 


To compute, for arbitrary x > 0 and arbitrary real a, 


the values of the incomplete gamma function, 


T(a,x) := f t e ie. (1) 
X 


2. Method 


Expansion in terms of reciprocal Laguerre polynomials. 


Let 


denote the Laguerre polynomial of order n with para- 
meter a (see ACCA I, p. 119). Then there holds for all 


x > 0 and alla < l 


T(a,x) = x e Dea a oe ee ee) 
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(see ACCA II, § 12.13). This expansion also holds for 
a 2 1, provided -x is not a zero of one of the polyno- 
mials ae see example [4] for how to deal with the 
exceptional case. Using the asymptotic theory of the 

Laguerre polynomials (see Szeg6, Orthogonal polynomi- 
als, Theorem 8.22.5) the n-th term of the above series 
(including the factor xe *) can be shown to be asym- 


ptotic to 
x _-4/nx 
an /% e 


The convergence of the series (2) thus is subgeometric. 
More precisely, the number n of terms required to get 


k correct decimal digits is asymptotically equal to 


2 
n= 0.33 a (3) 


pi 
and thus is roughly proportional to ae independent of 
a, and inversely proportional to x. 


For purposes of computation we write (2) in the 


form 
a ¥ 
(a,x) = Tr (4) 
n=l n-lon 
where 
Y _ Jo 5 x CE aa 
a A) eeu 


and 
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(5) 


From the recurrence relation of the Laguerre polynomi- 
als (see Szegd6, Orthogonal polynomials, equation 


(5.1.10)) we find 


hia = 0 , Loy =a 
L. = at (on -a-1+ x) -(n-a- 1) } 

n n n-l n-2 ° 
n= 1, 2, ... . TO save arithmetical operations, this 
is used in the form 

RL = Pita -a-—-1)(2 - & ) + (n + x)& }. (6) 
n n 1@ Weegee n-2 n= 1 
Summation is stopped if 
n 
<e€, 
nn-l 


where € may be chosen arbitrarily. 
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3. Flow Diagram 


Display 


T(a,x) := 8 


To save operations, most operations indicated in the 


n := n+ 1 box are anticipated in the box marked *. 
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4. Storage and Program 


Ro aT 2, R3 A a AG my 
"aeD ol : Yn le 
00 25. RCL 
01 CLX 26 RCL 1 
02 sfo 2 27 4 
03 sto 4 28 RCL 3 
04 1 29 $sTo 2 
05 sto 3 30 * 

06 STO 6 Sal + 
07 RCL 1 32 RCL 6 
08 RCL 0 33 : 

09 y* 34 + sTo 3 
10 sto 5 35 RCL 2 
11 RcL1 36 * 

12 e* 37 + 
13. sTOr5 38 «=6sTO+4 
+14 RCL 5 39 PAUSE 
15 RCL 6 40 ABS 
16 RCL O 41 RCL 7 
17 = 42 xX = y 
18 sTo*s 43 GTO 49 

19 1 44 1 

20 = 45 STO+6 

21 RCL 3 46 RCL 6 

o> REO 47 sTOr5 

23 7 48  GTo 14 


24 e + 49 RCL 4 
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5. Operating Instructions 


Load the program, move the operating switch to RUN. 
Select a mode of displaying results, for instance by 


pressing 


FIX 
9 
Load data as follows: 
a into Ro 
x i into Ry 
€ into Ro 
ee -10 
For 9 digit accuracy, let « = 10 . Pressing 
PRGM 
R/S 


will start computation. The calculator briefly dis- 
plays each term of the series (4) and stops if that 
term is < ¢€, displaying value of [(a,x). To get new 


value, it suffices to load new data and to press 


R/S 
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6. Examples and Timing 


The incomplete gamma function is of interest because 
it contains as special cases a number of important 

functions of mathematical physics and of statistics. 
The following special cases occur in the Handbook of 


Mathematical Functions by Abramowitz and Stegun. 


The exponential integral, 


20 at 
oe ap ane 
Xx 
Obviously is given by 
E, (x) = [(0,x) 
-10 
Using the tolerance e€ = 10 , the following va- 


lues are obtained for different values of x: 


number n computing 
of terms time 
x E, fx) required in sec 
5 0.559773594 67 190 
0 0.2193839 34 36 100 
2.0 0.048900511 19 54 
0 0.003779352 10 27 


The computed values are in error by at most one 


unit in the last place. 
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[2 | The exponential integrals of higher order, 


—-xXxt 
e 


n 
t 


E (x) := f 
1 


may be expressed as 


B(x). = es r(l-n,x). 


Some sample values computed by our program (using 


eE = 107”) are as follows: 

x E, (x) Eo E49 (%) 
0.5 0.3266439 0.1652428 0.0310612 
1.0 0.1484955 0.0860625 0.0183460 
2:0 0.0375343 0.0250228 0.0064143 
5.0 0.0009965 0.0007830 0.0002783 


All values are correct to the number of digits 


given. Computing times are comparable to those 
given under ; 


[3] The error function, 


> * 
eri (x), = Je adi 3 
yt O 


equals 
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as may be seen by a simple substitution. Using 


the tolerance ¢€ = 1 the program permits the 


following values to be obtained: 


number n computing 
of terms time 
x erf (x) used in sec 
1.0 0.842700793 34 97 
2.0 0.995322265 10 30 
3.0 0.999977910 5 15 


The values given are correct; however, the pro- 
gram Error function is much more efficient for 


Dy 


Chi-square probability function. This function 
may be defined as 


2 

(Ce es 

2 io? -9 
Q(x |v) := ‘ 
(3) 


Some values provided by our program (uSing € = 


10°) are as follows: 


Vv e 

1 2 3 4 
2 0.60653 0.36788 0.22313 0.13534 
4 0.90980 0.73578* 0.55783 0.40601 
6 0.98561 0.91970 0.80885 0.67669* 
8 0.99825 0.98101 0.93436 0.85712 
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As mentioned in section 2, for each a 2 1 the 
program may fail for isolated values of x due to 
negative zeros of the Laguerre polynomials. For 
instance, in the above table (x7 |v) = (2|4) is 
such a pair of values. Failure to compute [I (a,x) 
at x = Xo is recognized by an error halt. In such 
cases, [(a,Xo) may be obtained as the arithmetic 
mean of the values [(a,X_+9) where 6 is suffi- 
ciently small. In the above table, the values 
marked * were obtained in this manner. Of course, 
near failure may occur for values near Xo This 
is recognized by the appearance of large terms in 
the series (4). Again, more precise values can 


then be obtained by linear interpolation. 
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ERROR FUNCTION 


es Purpose 


To compute the error function, 
x 
fe dt , (1) 
0 


for all x 2 O with an error of less than nee 


2. Method 


A power series is used for 0 x < 2, and a continued 


> 
fraction for x = 2. 


a) The Taylor series x = 0, 
00 n 
=, 2 (-1) 2n+1 
erf(x) = 7 cae nAroontly * ; (2) 


is not to be recommended for numerical computation, be- 
cause the large alternating terms may cause cancella- 
tion, and because the general term is awkward to com- 


pute recursively. However, the series (2) may be ex- 


234 Special Functions 


pressed as a confluent hypergeometric series (ACCA I, 
GL) 
2 L 3 2 
erf(x) = —— X jF, (54 54 -*) - 
v1 
We thus may use Kummer's first identity (ACCA I, eq. 
(1.5-3); ACCA II, eq. (9.7-8)) to obtain 
2 oe 3 2 
eri (s)-= == 2 Fi (1 ; 5 # xX) . (3) 


In the last series, all terms are positive, and no 


cancellation can occur. Moreover, if 


2 00 
ex a i (13 : 7 xX ) = 2 ans 
v1 n=0 
then 
So. 2X a 
ron ee ’ 
: v1 
2 
a= * a 
n 1 “n-1 ” 
n + 5 
n= ly 2, This formula is used in the program to 


evaluate the series (3) for 0 = xX <2; 


b) Numerical experiments show that the series (3) 


can be used safely to compute erf(x) at least up to 
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x = 5.0, where erf(x) = 1 to the accuracy required. 
However, the time required to compute erf(x) by (3) is 
about 10x seconds for x 2 2 due to the increasing num- 
ber of terms that have to be summed. For x = 2 we 
therefore use an alternate method that is considerably 


faster. For all x > 0 there holds 
erf(x) = 1 - erfc(x) , 
where 


erfc(x) := ae 


VT 


is the complementary error function. This function ad- 


mits the asymptotic expansion 


Z 
-X 


erfc(x) % < 
v1 xX x 


which however cannot be used to obtain arbitrarily ac- 
curate values for any x. An expression that converges 
for all x > 0 is obtained by converting the asymptotic 
series into a continued fraction (ACCA II, § 12.13). 


This yields 


1 3 
-x L = | 5 | | 
erfc(x) = ip aid Fad Gad Gea (5) 


It has been determined experimentally that for x 2 2 
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the 16th approximant yields erfc(x) with an error < 


as Thus our program uses the approximation 


2 1 
erf(x) = 1 - os me + a ee sh 


Working with a fixed number of approximants has the 
advantage that the continued fraction can be evaluated 


rapidly working from the tail upward (ACCA II, § 12.1). 


3. Flow Diagram 


We set a :=a,S :=2Lt a,v?2=n + a m:= 8 — =k. 
n n 2 
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No 


No 

Display 
No 5 
erf(x) =1- a 


Display 


erf (x) 
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4. Storage and Program 


Ry Ry Ry R, Ry Fe Ry R, 
zg a Vv S m /T 5 E 

00 25  STo*l 
01 sto 0 26 RCL 1 
02 «7 27. =STO+3 
03 CHS 28° RCL 
04 e* 29 x<y 
05 RCL 5 30 cTo 18 
06 + 31 RCL 3 
07 sto l 32 GTO 00 
08 2 > 33. RCL O 
09 RCL O 34 8 
10 «zy + 35 sto 4 
11 GTO 33 36 lx Sy 
12. sTo*l 37 + 
13. RCL 6 38 RCL O 
14 +~sTo 2 39 + 
15 stot 40 RCL 4 
16 RCL 1 41 RCL 6 
ii 70'S 42 = 

+18 RCL O 43 x #0 
19 < 44 GTO 35 
20 RCL 2 45 R + 
21 1 46 sTOotl 
22 + 47 1 
23 sto 2 48 RCL 1 
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5. Operating Instructions 


Load the program; move the operating switch to RUN. 


Press 


FIX 


for nine-digit display (results will be accurate to 9 
digits). Load € := iors into Ro: 


Load 0.5 into R-: 


6 


STO 6 


Load /n into Re? 


VX 
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To obtain erf(x) for arbitrary x = 0, load x into 


X register and press 


PRGM 
R/S 


Value of erf(x) will be displayed after a few seconds. 


To obtain erf(x) for new x, load new x into X and press 


R/S 


OY 


. Examples and Timing 


1.99 yields erf(x) = 0.995111413, time 20 sec. 
= 2.00 yields erf(x) = 0.995322265, time 1llsec. 
0.01 yields erf(x) = 0.011283416, time 3 sec. 

5 ( dg 


000000000, time 11sec. 


~ Mm MR XM 
| 


-00 yields erf(x) = 


[el Lal [el [M1 [4 


The changeover point from power series to conti- 
nued fraction may be changed to any integer k 
such that 0 < k < 10, by substituting k for 2 in 
instruction 08. It is thus possible to compare 


the performance of the formulas (3) and (5) for 


various values of x. For instance, for x = 4 we 
get 
by (3) erf(x) = 0.999999985 , time 40 sec, 


by (5) erf (x) 0.999999985 , time ll sec. 
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COMPLETE ELLIPTIC INTEGRALS 


1. Purpose 


To compute, for 0 S k < 1, the complete elliptic inte- 


gral of the first kind, 


2. Method 


The arithmetic-geometric mean. Let Ob gs Bo be positive 


numbers, and let 
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rapidly converge to a common limit W(t 51 Bo)s called 
the arithmetic-geometric mean of Oy and By: 1 Be Ay = by 


Bo =k' := fl - e, then 


Tl 


rok) 


Moreover, if Yo f= ky 


(see ACCA II, § 9.10, problems 9 - 13). 


3. Flow Diagram 
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B := k' 
Oo 3= Ee 
aq 2= 1 
T ee 1 


Yes 
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4. Storage and Program 


k in xX 

00 25 RCL 7 
ot ean" 26 RCL 4 
02 cos 27 * 
03 STO 3 28 STO+6 
04 RCL 0 29 RCL 2 
05 a 30 RCL 3 
06 STO 6 31 * 
07 1 32 Vx 
08 STO 2 33 STO 3 
09 sto 7 34 RCL 5 

+ 10 RCL 2 35 STO 2 
11 RCL 3 36 GTO 10 
12 + 37 T 
13 2 38 2 
14 t 39 RCL 2 
15 sto 5 40 * 
16 RCL 3 41 + 
17 - 42 sT0 5 
18 < 43 1 
19 STO 4 44 RCL 6 
20 RCL 1 45 2 
21 xy 46 : 
22 GTO 37 47 - 
23 2 48 * 
24 STO*7 49 RCL 5 
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5. Operating Instructions 


Load the program; move the operating switch to RUN. 
Select a mode of displaying numbers, for instance by 


pressing 
FIX 9 


to get 9 digits after decimal point. Load af = 10. 


into R as follows: 


1 3f 


STO l 


Load k into Ro and also leave it in X register. Pres- 


sing 


PRGM 
R/S 


will after a few seconds cause K(k) to be displayed. 


To display E(k), press 


m* 
AV 
Ke 


K(k) then will be in Y register. 
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6. Examples and Timing 


For k = 1/V2 = 0.707106781 we get 


K = 1.854074677 , E = 1.350643881 ; 


computing time approx. 8 sec. For k = 0.9999 the re- 


sults are 


K = 5.645148167 , E = 1.000514502 ; 


computing time approx. 12 sec. 
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BESSEL FUNCTIONS, INTEGER ORDER 


1. Purpose 


To compute the Bessel functions J, 6) for arbitrary 


x > 0 and for the orders n= 0, 1, 2, 3, 4. 


2. Method 


J.C.P. Miller's method of backward recurrence [see W. 
Gautschi, Computational aspects of three-term recur- 
rence relations, SIAM Review 9, 24 - 82 (1967)]|. We 
use a recurrence relation satisfied by the sequence 


{J (x) J (x fixed) in the form 


_ 2n 
eae es a ers (1) 
If the recurrence 
2n 
pat a” a 2 
real x Yn Ya (2) 


is solved with arbitrary starting values, say, 
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Ym+1l *~ r Ym tT 1, 


where the integer m is sufficiently large, then the 
elements of the sequence y soon become 


J 


m-n’ ~m-2’ 


proportional to J .. , because {33 isa 


m-1’ ~“m-2’ 
"recessive" solution of the difference equation (2). 


Thus, approximately, 


1 
J eae: ’ (3) 


where c is independent of n. The factor c may be de- 
termined from the known relation 


Jy (x) + 23, (x) + 273,(x) + ... = 1, 


4 
which gives 


c= Yo + 2¥5 + 2Y4 in ee (4) 


In this program m may be selected as an arbitrary 
positive integer; see the examples for a rule of thumb 
on how to select m as a function of x. The program 
then constructs the sequence Lye always saving the 
last five elements. The values Js) are then computed 


from (3), approximating c by 


Mg, EN ge Ng See POLY rao ° 
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3. Flow Diagram 
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4. Storage and Program 


NO FPF WwW NO FSF W © 
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5. Operating Instructions 


Load the program; turn the operating switch to RUN. 
Select the mode of displaying numbers, for example, 


by pressing 


SCI 8 
Load m into R.. (m = 30 is adequate for values x < 10.) 
Load x into X. Press 
PRGM 
R/S 


to start the computation. The calculator will stop by 


displaying calculated value of J,(x); the values J, (x) 


0 


(k = 0, 1, ... , 4) are also found in Ry To restart 
the computation with a new value of x, simply reload 


m into R the new x into X, and press 


7 ¢ 


R/S 


6. Examples and Timing 


For x = 1, m = 30 we get the values 


0.76519769 
0.44005059 


GQ 
oO 
II 
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J. = 0.11490348 
J, = 0.01956335 
Jy = 0.00247664 


which are exact to the number of digits given. 


Computing time 40 sec. 
For x = 10, m = 30 the results are 


= -0.24593576 
= 0.04347275 
= 0.25463031 
= 0.05837938 
-0.21960269 


correct to all digits given. 


For x = 17.5 we choose m = 40. The computing time 


now is about 59 sec. There results 


Jo = -0.10311040 
Jy = -0.16341997 
J. = 0.08443383 
J, = 0.18271913 
J, = -0.02178727 
Again, the errors are < io. The value m = 50, 


which theoretically produces still more accurate 


results, causes overflow. 
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BESSEL FUNCTIONS, ARBITRARY ORDER 


1. Purpose 


To compute, with an error < 16>", the values of the 


Bessel function J. (x) of arbitrary real order v for 
< < 
O =x = 10. 


2. Method 


Power series expansion. We have (see ACCA II, § 9.7) 


2 
x, V n,xX nN 
oe ue : eis c 
Vv Vv n=0 Vv wot 


We write this in the forma } bo where 
v 
( 


x 
5) 
rt b = AL 


ais T(v4l) ° 0 ’ 


and the remaining b are calculated recursively by 
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The summation is terminated if |b. < 10°**. the term 


[T(v+l) required in the computation of a must be compu- 
ted separately, if necessary by the program gamma 


function. 
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3. Flow Diagram 


We write s := 2} b 


No 
ib! < 10 ? n := n+l 


Yes 


Display J.) (x) = ab 


256 


4. 
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Storage and Program 


0 Ry 9 

Peery | ao. * 
00 
01 2 
02 + 
03 sto 7 
04 1 
05 EEX 
06 CHS 
07 it 
08 2 
09 SsT0O 2 
10°. RCL? 
11 RCL O 
12 y* 
13 RCL 1 
14 2 
15 SsT0o 3 
16 1 
17 sto 4 
18 sTo 5 
19 STO 6 

+ 20 RCL 7 
21 eP 
22 RCL 6 
23. RCL O 
24 + 


NOlLK 
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5. Operating Instructions 

Load the program, and move the switch to RUN. Press 
FIX 8 

to get a reliable eight-digit output. Load 


5 : 
into Ro 


r (v+_1) into Ry 
(The latter value may be computed by the program gamma 


function.) Load x into the X register. When 


PRGM 
R/S 


is pressed, the calculation will start. The calculator 
stops by displaying J (x). 
This program does not work if v is a negative in- 


teger, v = -n. Here the relation 
n 
Ti Ax) = AGL) de) 


should be used. The program may be inaccurate if v is 
close to a negative integer. 
For nonintegral v, the Bessel function of the se- 


cond kind may be computed using the relation 


258 


For integer values of 
going aS v > n may be 


the Romberg algorithm 


6. Examples and Timing 
J) (2) = 0.2238907 
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(x) cosvt - J_ (x) 

sinvtT 
vv = n, the limit of the fore- 
calculated numerically, using 
(see example [4]). 
8, all digits correct, time re- 


quired about 13 sec. 


J,(10) = -0.24593576, all digits correct, time 
required about 27 sec. 
To compute J379(5) , we require 
5. 1 3 
rig) = 5° 5 v0 
Computation yields J372©) = -0.16965131, in 


agreement with the exact value. 
To compute Yo (x), we may use the relation 


2 oO 
T OV Jes) v=0 


Letting 
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our program yields for x = 2 
Vv qd. 
2 0.00000000 
1 0.57672481 
Ces: 0.74780184 
0525 0.78844828 
0125 0.79839964 


With the Romberg algorithm the limit of these va- 
lues as v > 0 is determined as 0.80169622. Thus 


we obtain 
Y,(2) = 0.51037566 , 


which is in error by only l * ie. 
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BESSEL FUNCTIONS: ASYMPTOTIC SERIES 


de Purpose 


To compute the asymptotic series, for arbitrary x > 0, 
of the Bessel functions of both the first and second 


kind, 


for arbitrary real orders v. 


2. Method 


The asymptotic series is computed most conveniently 
from the formula 
page 


i(x ) 
Ain f D 2 4 l l 1 
J (x) +i Y,, &) /Z e Fo (5 + V5 - Vv; Sin)! (1) 


where is a hypergeometric series (see ACCA II, 


2F 0 
§ 11.5). We write 
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and calculate the a by means of the recurrence 


1 1 
(- 5 + v + n)(- 5 7 CG + n) 


0 n : 2xn n-l 


The sums R and S are then given by 


The point of termination of these series may be chosen 


arbitrarily. Ordinarily, to get best results, one 


should terminate as soon as }a,| > |a . If v equals 


an integer plus Y2, the sums R and S aa terminate 
mathematically, and the representation given by (l) is 
not only asymptotic but exact. 

To save programming storage, the correct assign- 
ment of the terms ta to the sums R and S is left to 
the operator (see operating instructions). It then be- 
comes possible to compress the whole program into 49 


steps by converting R + iS into polar coordinates, 


R+i1S = Te (ES OO} sy 


which yields 
VIO TT VOT TT 
re ae ay een De ag ee 
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Converting back to carteSian coordinates then yields 


J (x) ay [= T cos(x - > - 7 +O) 5 (2) 
¥ (x) av . T sin(x - a Soe Beis) (3) 


3. Flow Diagram 


We set q := (-S4+vtn) (-5-vtn) /2xn. 
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n+1 


Display q, a 


Summation to 


be continued? These deci- 


sions to 


be made 


Yes 


manually 
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4. Storage and Program 


PAUSE 


AVI 


Bessel Functions: Asymptotic Series 265 


5. Operating Instructions 


Load the program, and switch to RUN. Press 
FIX 8 


to obtain adequate display of results. Press 


to cause arguments of trigonometric functions to be 


interpreted in radians. Load 


Xx int 
1 O Ro 


2v+l into Ry 


(the index v must be stored in this form for reasons 


Of space). Pressing 


PRGM 
R/S 


will first cause 2 to be computed and to be displayed 
briefly. If Iq, | < 1, allow computation to continue by 
doing nothing. The machine will stop by displaying a. 


Press 


4k+1 
4k+2 


STO-5 if n 
STO-4 if n 
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STO+5 if n 4k+3 


STO+4 if n 


4k+4 
Pressing 
R/S 


will cause the computer to proceed and to calculate 


next qn and next ane If Iq, 2 1, quickly press 


R/S 
GTO 29 
R/S 


At next stop, the calculator will display asymptotic 
value of J (x). To get asymptotic value of Y Ox), 


press 


AN 


Caution: If R < 0 and |S/R| < toma the value of 
will be in error by precisely 1, and signs of both J, 


and x will be wrong. 


6. Examples and Timing 
[1] v = 0, x = 10. Stopping at n = 20 yields 


J, (10) % -0.24593576 , YQ (10) % 0.05567117 
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Both values agree with the exact values to all 


digits. Total computing time about 95 sec. 


J,,(2) = 0.51301614 , Yo (4) = 0.23478571 


These are the exact values. 


) = 0.19056437 , ) = =-0.57174942. 


yo? Se aa 


These values agree with those computed by the 


program Bessel functions, arbitrary order. 


268 Special Functions 


RIEMANN ZETA FUNCTION ON CRITICAL LINE 


1. Purpose 


To compute the function 


E(t) := nS + it) 


for real values of t, where 


ECs). 2] 8 m (Re s > 1) 


The function € is entire, even, and real for t real. 


There holds €(t) = 0 if and only if s = - + it isa 


nontrivial zero of t(s) (see ACCA II, § 10.8). 


2. Method 


We use the integral representation 
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E(t) = sf £(x) cos(2tx) dx , (1) 
0 
where f(x) := o"(x) - (x), 
d(x) i= e* mea ’ 
Kay ae —n Tu 
n=1 


(see Bieberbach, Funktionentheorie Vol. II, Chapter 8). 
The integral (1) is approximated by the midpoint 


formula, uSing an arbitrarily chosen step h, that is, 


(<0 ) 
E(t) =h x £(x,) cos(2tx,) , (2) 
k k 
k=0 
where xy := sh + kh. The summation is terminated as 
soon as J £(x,) | < ¢€. The function f is evaluated by 
means of the series 
(00 ) 7 4x 
See oS Weare Once ee (3) 


The summation is stopped when the general term (which 
is always positive) becomes < ¢€«. The series converges 
rapidly for all x 2 O. 

Because Of lack of programming space, no automatic 


step refinement is built into the program. However, 


270 Special Functions 


the step should always be refined manually to check 
the accuracy of values. Because f(x) = O(exp ga") 
as x > ~, the error of the numerical values of the 
integral (1) is O(h") for every m > 0, and Romberg 
acceleration is therefore neither required nor effec- 


tive. 


3. Flow Diagram 


We denote the sum in (2) by &, the sum in (3) by S, 
and the n-th term of S by ane 
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z 
2 
0 


2= n+l 


S 


Yes 


Display €(t) := mr s= 2 + £(x)cos(2tx) 
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4. Storage and Program 


RCL 3 
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oe Operating Instructions 


Load the program, and switch to RUN. Press 


SCI 8 


to get the floating eight-digit display of results. 
Load the data as follows: 


2t into Ro 
h into Ry 
€ into R 
l 2 
5h into R ‘ 
Clear Ry by 
CLX 
STO 4 


Because argument of cosine is in radians, press 


When 


PRGM 
R/S 


is pressed, the calculator will start whirling and 
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stop by displaying &(t) as computed by formula (2). To 


check accuracy, the computation should be repeated 


with h := sh. The necessary new data are furnished by 
RCL 1 2 
2 + 
= STO 3 
STO 1 CLX 
STO 4 


Note: Some of the foregoing instructions can be 
automated on a calculator with a few more memory lo- 


cations. 


6. Examples and Timing 


=L2 


Using ¢« := 10 , we get for 

t = 14 h= 0.2 — € = -0.4648/7914 ( 24 sec) 
Oe: +0 .00063699 ( 49 sec) 
0.2.05 0.00020129 ( 99 sec) 
0.025 0.00020129 (198 sec) 

t = 15 h= 0.1 —~ € = -0.00001633 
0.05 -0.00070570 
O3:025 -0.00070570 


Thus a zero exists between t = 14 and t = 15, implying 
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the existence of a zero of the Riemann zeta function 


on the critical line between s = : + 141 and s = ; + 


2 
15i. 
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Aitken acceleration, 48, 53 
asymptotic expansion, 162, 235 


for Bessel function, 260 


Bell number, 150 
Bernoulli number, 132 
Bernoulli's method, 80, 87 
Bessel function, 
of arbitrary order, 253, 260 
of integer order, 201, 247 
of second kind, 257, 260 
binary expansion, 126 


binomial coefficient, 18, 72 


chi-square probability function, 23l 
continued fraction, 24, 31, 235 


convergence, quadratic, 53 


deflation, 74, 77 
difference, finite, 196 
differential equation, 
autonomous, of second order, 189 
of first order, 181 
linear, of second order, 189 
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278 Index 


elliptic integral, complete, 241 
Enestré6ém-Kakeya theorem, 118 

equation, 41 

error function, 230, 233 

Euclidean algorithm, 14, 19, 33, 34, 36 
Euler's constant, 167, 175 

exponential integral, 229 


of order n, 230 


fixed point, 41, 53 


Fourier transform, fast, 124 


gamma function, 214 
incomplete, 223 
Gauss-Lucas theorem, 113 


greatest common divisor, 17, 19 


Horner algorithm, 67, 73 
group property of, 72 


integration, numerical, 155, 166 
irrationality, quadratic, 30, 3l 


iteration, 4l, 47, 53, 182 
Lagrange-Blirmann theorem, 140 
Laguerre polynomial, 223 


Log-Arcsine algorithm, 207 


map, contracting, 42 
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mean, arithmetic-geometric, 241 
midpoint formula, 155, 269 
Miller, J.C.P., formula of, 135 


Newton's method, 
for complex roots, 59 
for polynomials, 73 


for reciprocal power series, 124 


pendulum, 186, 194 
Plana summation formula, 170 
polynomial, stable, 101 

with zeros in unit disk, 107, 114 
power series, 

exponentiation of, 144 

power of, 134 

reciprocal of, 123 

reversion of, 140 


prime factor, 9 


quotient-difference algorithm, 


87, 95, 216 


recurrence, backward, 247 

Riemann zeta function, 177 
on critical line, 268 
zeros of, 275 

Romberg algorithm, 156, 160, 161, 162, 169, 
LPL: LG, LIL y. LTBy 1867-183, 195. 209.259 
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Routh algorithm, 101 
Runge-Kutta method, 182, 189 


Schur-Cohn algorithm, 107, 114 
Stirling's formula, 142, 214 


trapezoidal rule, 181 


